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Abstract 

Micro-propulsion  has  been  studied  for  many  years  due  to  its  applications  in  small-to- 
medium  sized  spacecraft  for  precise  satellite  attitude  control.  Electrospray  thrusters 
are  promising  thrusters  built  upon  the  state  of  the  art  in  micro-technology  and  with 
flexible  performance  in  terms  of  their  high  efficiency  and  high  specific  impulse.  One 
challenge  is  to  investigate  in  detail  the  mechanism  for  ion  emission  to  complement 
experimental  results  and  understand  better  how  emission  occurs  in  the  micro  to  nano 
scale.  Thus,  atomistic  modeling  is  used  to  understand  properties  of  emitted  charged 
particles  which  determine  how  the  thrusters  perform. 

As  a  preliminary  study  of  ion  emission  from  Taylor  cones,  ion  evaporation  from  3 
-  5  nm  droplets  was  observed  in  molecular  dynamics  (MD)  simulations  to  validate  the 
atomistic  modeling  and  to  investigate  activation  energies.  Ion  emission  was  examined 
in  terms  of  internal  and  external  electric  fields  and  the  activation  energies  of  each  case 
were  obtained  using  Schottky’s  model  and  direct  energy  calculation  to  compare  with 
experimental  values.  Ion  emission  was  mainly  observed  with  electric  field  strengths 
between  1.2  -2.0  V/nm  and  the  emitted  species  include  both  solvated  and  non-solvated 
ions.  Propulsive  properties  from  Taylor  cones  are  examined  using  results  from  the 
analysis  of  electric  current  from  ion  emission. 

In  addition  to  an  observation  of  ion  emission  from  liquid  droplets,  numerical  simu¬ 
lations  for  interactions  between  a  solid  plate  and  liquid  droplets  were  conducted  with 
MD  simulation.  It  was  concluded  that  another  selection  of  force  field  needs  to  be 
considered  to  pursue  further  details,  such  as  electrochemical  effects. 

Thesis  Supervisor:  Prof.  Paulo  C.  Lozano 
Title:  H.  N.  Slater  Assistant  Professor 
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Chapter  1 


Introduction 

1.1  Electrosprays 

Electrospray  is  a  technology  to  produce  ions/droplets  from  conductive  liquids  using 
a  relatively  simple  setup  as  shown  in  Figure  1-1.  Ion/droplet  emission  occurs  at  the 
apex  of  a  liquid  cone  which  is  generated  when  a  high  electric  potential  between  an 
extractor  and  conductive  liquid  is  applied.  There  are  several  types  of  emitters  with 
the  most  common  being  the  capillary-type  often  used  in  mass  spectrometry.  The 
other  types  are  externally  wetted  emitters  (Figure  1-2)  and  porous  emitters.  For  any 
kind  of  emitter,  there  are  two  underlining  requirements:  to  locally  enhance  electric 
fields  at  the  tip  so  that  a  liquid  cone  appears  and  to  transport  the  liquid  from  a 
resevoir  to  the  emission  site.  The  cone  is  called  a  Tayor  cone  in  honor  of  G.I.  Taylor 
who  first  derived  a  mathematical  model  for  ideal  conditions  in  1964  [1],  Details  of 
the  emission  mechanism  and  Taylor  cones  will  be  discussed  in  section  3.1  and  3.2. 

As  for  the  applications  of  elecrosprays,  they  can  be  used  in  a  broad  field  of  studies 
from  biology  to  space  technologies.  Mass  spectrometry  for  biological  molecules  is 
the  most  well  know  application  of  electrosprays  and  was  subject  of  the  Nobel  prize 
winning  work  by  J.Fenn  in  2002  [2],  Some  other  technologies  are:  thin-film  deposition 
and  painting  [3],  electrospining  which  is  used  to  create  uniform  strands  of  nano-scale 
fibers  [4],  micro  scale  atomizers  for  fuel  injection  [5], [6],  molecular  imaging  applica¬ 
tions  [7],  and  focused  ion  beams  (FIB)  [8], [9], [10]. 
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Figure  1-1:  Schematic  of  electrospray  setup 


Figure  1-2:  Tip  of  externally  wetted  electrospray  thruster.  Courtesy  of  A.  Zorzos 
(MIT) 
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An  application  for  space  propulsion  is  explained  in  chapter  2. 


1.2  Electrospray  History 

The  use  of  electric  fields  for  controlling  fluids  can  be  traced  to  the  eighteenth  and 
nineteenth  century.  [11], [12]  The  breakup  of  a  charged  liquid  drop,  Coulomb  fission, 
was  predicted  by  Rayleigh  in  his  analytical  work  in  1882.  [13]  He  presented  that  the 
instability  occurs  when  the  forces  of  electrostatic  repulsion  of  surface  charges  exceed 
the  forces  of  surface  tension  for  a  droplet.  Zeleny  first  conducted  an  experiment  of  the 
observation  of  the  Rayleigh  instability  and  visual  formation  of  cone-jet  structures  with 
a  diluted  solution  of  hydrochloric  acid  in  1914.  [14], [15]  However,  the  research  was 
hidden  from  space  propulsion  applications  for  a  half  century  until  Krohn’s  electrospray 
research  in  1961.  [16], [17] 

Unfortunately  electrosprays  were  not  pursued  further  as  a  thruster  in  the  1960’s 
because  the  value  of  mass  to  charge  ratio  (m/q)  was  too  large  which  requires  high 
accelerating  voltages  (10  to  lOOkV)  from  V  =  rj^~,  even  though  higher  particle  mass 
is  an  attractive  feature  because  it  can  effectively  increase  the  thrust  density.  [18]  In 
the  same  period,  the  theoretical  model  of  cone  formation  was  developed  by  Taylor 
(see  detail  in  section  3.2).  [19]  Since  then,  electrospray  research  did  not  see  significant 
progress  until  the  breakthrough  achieved  by  Fenn  in  the  late  1980’s,  escalating  an 
interest  in  nanotechnology.  Fenn’s  contribution  made  it  possible  to  utilize  electrospray 
technology  in  mass  spectrometry,  while  at  the  same  time  enabling  a  reduction  of  the 
mass  to  charge  ratio.  Fenn  was  awarded  the  Nobel  Prize  in  Chemistry  in  2002.  [20] 
With  this  achievement  and  growing  interest  in  ionic  liquids,  electrospray  propulsion 
research  embarked  on  a  new  era.  Ionic  liquids  are  molten  salts  at  room  temperature 
containing  anion  and  cations  exclusively,  which  allow  a  purely  ionic  regime  of  emission 
during  electrospray  use  for  extraction  voltages  under  2k V.  [21], [22], [23] 

As  for  studies  of  ion  evaporation  from  a  droplet,  Thomson  and  Iribarne  derived 
the  threshold  electric  field  at  which  ion  field  evaporation  occurs  before  Coulomb 
fission.  Emission  is  triggered  by  adding  the  energy  to  overcome  the  barrier  associated 
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with  both  the  solvation  energy  of  the  evaporated  ion  and  the  electrostatic  interaction 
between  the  ion  and  the  droplet  with  the  assumption  that  the  barrier  is  approximately 
4-5  A  from  the  droplet  surface  [24].  In  addition,  they  showed  that  the  critical  electric 
field  is  approximately  1  V/nm  for  droplets  less  than  10  nm  in  diameter  using  a  kinetic 
model  and  charge  mobility  experiments  in  the  late  1970’s.  [25]  The  Loscertales  and 
Fernandez  de  la  Mora  group  and  the  Fenn  group  concluded  the  critical  value  for  the 
evaporation  of  singly  charged  ions  was  between  0.7  and  1.9  V/nm.  [26] 

The  detailed  mechanisms  of  field  evaporation  are  still  not  well  understood  because 
the  complex  particle-particle  interactions  of  the  thermalized  droplets. 
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Chapter  2 


Electrosprays  in  Space  Propulsion 


2.1  Electric  Propulsion 

Space  propulsion  systems,  and  rockets  in  general,  provide  force  to  accelerate  space¬ 
craft  consuming  propellant.  [27]  The  propulsion  systems  are  used  both  for  launch 
of  spaceship  from  the  Earth  and  for  moving  spacecraft  in  space.  For  instance,  space 
thrusters  are  used  for  orbital  maneuvers,  north-south  station  keeping  (NSSK),  orbital 
altitute  and  attitude  maintenance  and  interplanetary  travel  for  deep  space  explo¬ 
ration.  There  are  mainly  two  types  of  the  systems:  chemical  and  electric  (Figure. 2-1). 

Chemical  propulsion  thrusters  work  by  the  conversion  of  an  enthalpy  increase 
from  the  energy  of  chemical  reactions  into  jet  kinetic  energy  of  reaction  producs 
(combustion),  while  electric  propulsion  makes  use  of  electric  energy  to  produce  thrust. 


Chemical 

Electrical 


Augmented 

Electrostatic 

Electromagnetic 

Radio 

Frequency 


> 


L-> 


Electrospray 

Hall  Thruster 
Ion  Thruster 


Figure  2-1:  In-space  propulsion  technology 
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Electric  propulsion  also  needs  propellant  to  accelerate  and  since  thrust  depends  on 
the  electric  power  available,  the  thrust  of  electric  propulsion  is  typically  smaller  than 
chemical  propulsion.  Electric  propulsion  is  used  in  space,  a  vacuum  environment,  due 
to  the  low  densities  required  by  most  electric  discharges.  [28] 

There  are  several  ways  to  transfer  electric  power  to  the  propellant:  electrothermal, 
electrostatic,  and  electromagnetic  devices. 

Electrothermal  thrusters  include  the  Resisto-jet  and  arc-jet  which  are  augmen- 
tated  devices  which  are  primarily  chemical  thrusters  but  use  electrical  power  to 
improve  their  performance.  The  propellant  gas  is  heated  electrically  and  expanted 
through  a  nozzle  to  convert  its  thermal  energy  to  a  jet  of  directed  kinetic  energy.  Elec¬ 
trospray  thrusters  are  electrostatic  devices  as  well  as  ion  thrusters  in  which  the  electro¬ 
static  forces  accelerate  ions/droplets.  [29]  Working  details  of  electrospray  thrusters 
are  explained  in  section  2.1.1.  For  electromagnetic  thrusters,  the  acceleration  of  a 
body  of  ionized  gas  is  obtained  by  the  interaction  of  currents  driven  through  the  gas 
with  magnetic  fields  established  either  by  the  currents  or  by  external  means. 

2.1.1  Electrostatic  Thrusters 

Typical  electrostatic  thrusters  are  ion  thrusters  and  electrospray  thrusters.  In 
these  devices,  the  electrostatic  forces  that  accelerate  the  ions/droplets  are  also  directly 
applied  to  electrodes,  and  this  is  how  the  structure  receives  thrust.  Ion  thrusters 
extract  ions  from  ionized  propellant  by  a  potential  drop  through  grids  while  electrons 
are  emitted  from  a  cathode  which  is  placed  at  the  exit  of  the  thruster.  These  electrons 
neutralize  the  ion  beams,  avoiding  spacecraft  charging.  Ion  thrusters  typically  have 
specific  impulses  in  the  3000  to  4000  range  with  power  efficiencies  as  high  as  75 
percent.  One  of  restrictions  in  ion  thrusters  is  space  charge  limitation  which  is  a 
phenomena  in  which  the  amount  of  extracted  ions  are  limited  due  to  their  shield  effect 
to  the  applied  fields.  On  the  other  hand,  electrosprays  are  not  largely  affected  by  space 
charge  limitation  because  the  emitter  to  extractor  spacings  are  significantly  smaller  in 
electrospray  devices.  In  addition,  electrospray  thrusters  do  not  suffer  from  problems 


22 


associated  to  gaseous  propellants.  One  of  the  main  issues  in  electrospray  thrusters 
is  the  transportation  of  liquid  propellant  from  a  reservoir  to  the  emitter  needle  tips. 
In  recent  work,  micro-fabricated  porous  emitters  have  been  studied  as  they  provide 
sufficient  hydraulic  impedance  to  sustain  emission  even  at  low  currents.  [30], [31] 


2.2  Electrospray  Thruster  Application 

In  a  broad  sense,  field  emission  electric  propulsion  (FEEP)  and  colloidal  engines  have 
identical  principle  of  operation  with  our  electrospray  thrusters.  The  main  differences 
are:  (a)  FEEP  uses  a  liquid  metal  which  requires  high  operation  voltages  (>5kV) 
due  to  its  high  surface  tension  and  produce  only  positive  beams.  [32]  Colloidal  engines 
can  be  used  with  organic  liquids  with  relatively  high  conductivity  and  lower  surface 
tension  which  allows  lower  voltages  (lkV).  (b)  Electrospray  and  FEEP  thruster  focus 
on  ion  emission  to  obtain  high  efficiency  although  colloidal  engines  could  be  nearly 
as  efficient  operating  with  droplets.  In  short,  electrospray  thruster  families  in  general 
can  adapt  to  missions  through  several  regimes  and  types  of  liquids.  Droplet  regime 
gives  low  Isp/current  due  to  its  large  specific  charge,  on  the  other  hand,  ion  regime 
gives  high  Isp/current  though  the  mixed  regime  decrease  thrust  performance.  [33] 
These  flexibility  of  characteristics  within  electrospray  thrusters  give  options  to  many 
kind  of  space  missions. 

The  following  is  a  list  of  potential  applications  for  electrospray  thrusters: 

1.  Precise  attitude  control  of  spacecraft 

The  small  thrust  allows  accurate  positioning  of  spacecraft.  One  example  is 
the  LISA  Pathfinder  mission  which  features  both  FEEP  and  colloidal  thrusters 
and  plans  to  be  launch  in  2013.  [34]  The  colloidal  thruster  was  developed  by 
Busek  Company  in  the  U.S..  One  of  advantages  of  electrospray  thrusters  is  that 
they  do  not  need  conditioning  time  which  is  usually  required  for  other  electric 
propulsion  and  FEEP  thrusters.  Liquid  metal  used  for  FEEP  propellant  needs 
to  be  heated  up  to  liquefy. 


23 


2.  Main  propulsion  for  Nano/Micro  satellites 

With  growing  demand  of  miniaturization  of  spacecraft  components  and  satellites 
themselves,  electrospray  thrusters  have  the  potential  of  providing  main  propul¬ 
sion  to  small  satellites.  CubeSats  are  one  example  of  small  satellites,  which 
are  widely  considered  in  universities  and  companies  because  of  their  advantages 
in  small  satellite  space  science  and  exploration  at  lower  costs.  CubeSats  and 
other  nano  satellites  have  typically  a  volume  of  1  liter  and  weight  no  more  than 
1  kg.  Most  thruster  technologies  are  larger  than  the  satellite  itself.  Also,  the 
fuel  is  fed  by  capillarity  action  in  electrospray  thrusters  therefore  avoiding  a 
complex  fuel  feed  system  which  takes  space  in  satellites.  Another  option  might 
be  micro-scale  cold  gas  thrusters  [35],  however  it  has  a  problem  to  miniaturize 
the  thruster  without  downgrade  the  performance  and  the  Isp  is  limited.  FEEP 
and  colloid  thrusters  have  a  possibility  to  contaminate  a  spacecraft  and  specific 
impulse  of  colloid  thrusters  is  generally  lower  than  electrospray  thruster. 
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Chapter  3 


Electrospray  Physics 


3.1  Surface  charge 


The  basic  physics  of  electrospray  can  be  expressed  by  balance  between  surface  traction 
by  external  electric  field  and  surface  tension  of  liquids.  The  external  electric  field  can 
be  related  with  the  electric  field  in  the  liquid  using  Gauss’  Law  [29],  and  for  dielectric 
liquid, 


En, liquid  —  E n,gas 


(3.1) 


I 

^n,  gas 

£0.  e  =  1 

gas 

e 

dee 

^  £0£ 

T  liquid 

I  ^n,  liquid 

Figure  3-1:  Schematic  of  surface  charge  for  dielectric  liquid 


Here,  e  is  dielectric  constant  of  liquid.  This  can  take  large  values  and  for  perfect 


25 


conductors,  the  field  inside  liquid  becomes  zero  which  is  caused  by  high  mobility  free 
electrons  with  fast  relaxation  times. 

Thus,  the  relaxation  time  can  be  considered  as  another  important  aspect  of  elec¬ 
trosprays.  The  liquid  takes  time  to  respond  to  external  electric  fields  and  this  depends 
on  the  mobility  of  charged  particles. 

Conductivity  is  expressed  using  mobility  of  positive  and  negative  particle: 

k  =  ne(/r+  +  n~)  (3.2) 

The  rate  of  free  surface  charge  density  (07)  accumulation  is, 

^  liquid 

where  the  surface  charge  density  can  be  obtained  by  Gauss’  law: 


(3.3) 


&  f  —  b)  En^gas  eco  En:iiqUid 

Prom  equation  3.3  and  3.4, 


(3.4) 


do/ 

dt 


—af  =  -En 

66q  6 


Solving  the  differential  equation, 


(3.5) 


a,  = 


(3.6) 


where  r  =  ^  is  the  relaxation  time,  or  the  time  to  reach  equilibrium  in  the  charge 
distribution. 


3.2  Taylor  cone 

Taylor  cones  are  produced  when  conductive  liquids  are  placed  under  an  external 
electric  field.  [29], [19]  Taylor  estimated  the  shape  assuming  that  the  cone  formed  by 
a  balance  between  electric  pressure  and  surface  tension.  This  is  expressed  by: 
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Figure  3-2:  Taylor  Cone  Schematic 


2C0 El, liquid  ~  (3-7) 

Here,  7  and  1  / Rc  represent  surface  tension  and  curvature  which  is  shown  in  fig¬ 
ure  3-2.  Using  geometric  relations,  the  right  hand  side  of  equation  3.7  is 

—  -  ~-cosa  =  —cota  (3.8) 

Rc  R  r  v  ' 

From  both  equations  3.7  and  3.8,  the  electric  field  along  the  surface  of  the  liquid 

cone  is  obtained: 


-'ll, liquid 


27 cota 
e0r 


(3.9) 


The  schematic  of  electric  field  along  liquid  cone  surface  is  shown  in  figure  3-3. 
From  equation  3.9,  it  is  obvious  that  electric  field  is  proportional  to  r~1/2  therefore  it 
is  smaller  upstream  and  at  the  tip  of  cone  the  electric  field  becomes  infinite.  Large 
electric  fields  induce  ion  emission.  However,  it  is  challenging  to  numerically  analyze 
the  tip  region  due  to  its  small  dimension  and  dramatic  potential  gradients,  thus 
alternate  models  are  used,  such  as  the  sphere  on  cone  (SOC)  model.  The  detailed 
calculation  is  shown  in  section  7.3. 

Following  his  calculation,  Taylor  mathematically  derived  that  the  inner  half-angle 
is  independent  of  liquid  properties  and  has  a  universal  value  of  49.3°.  Real  experi- 
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Figure  3-3:  Schematic  of  liquid  cone  and  electric  of  liquid  surface 

mental  cones  show  a  remarkable  resemblance  to  Taylor  cones,  even  though  the  apex 
singularity  in  real  cones  is  removed  by  the  emission  of  a  charged  liquid  jet  in  most 
solvents  [36]  or  the  emission  of  ions  in  ionic  liquids  and  metals. 


3.3  Activation  energy 

Ion  emission  can  be  explained  using  arguments  from  statistical  mechanics  and  defining 
the  energy  needed  for  ion  evaporation,  or  activation  energy.  Iribarne  and  Thomson 
first  presented  a  model  to  describe  ion  evaporation  from  charged  droplets.  [24] 

Here  we  derive  the  fundamental  equation,  Schottky’s  equation,  using  the  equilib¬ 
rium  between  ions  in  the  liquid  and  gas  phases.  [37]  In  equilibrium,  ions  evaporate 
from  the  liquid  surface  into  the  gas  phase  and  also  the  opposite  happens.  Equilibrium 
can  be  expressed  as: 


(/h)i,s  —  i)g,v  (3.10) 

where  /q  i,  l,  s,  g,  v  correspond  to  chemical  potential,  ion,  liquid,  surface,  gas  and 
volume. 

Assuming  Maxwell-Boltzmann  gas,  the  above  equation  would  be, 
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kTiln 


(3.11) 
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where  k  is  Boltzmann  constant  and  T  is  temperature,  which  in  equilibrium,  Tg  = 
Tt.  Q  and  N  represent  the  partition  function  and  number  of  molecules,  respectively. 

Considering  dimensional  aspects  which  Q  and  N  have  (e.g.  2D  for  liquid  surface 
and  3D  for  gas  phase),  the  equation  3.11  is, 
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(3.12) 


(3.13) 


Here,  A  is  area  of  the  liquid  surface  and  V  is  volume  of  the  gas.  Thus,  q  and  n  are 
the  partition  function  and  number  of  molecules  per  unit  area  and  volume,  respectively. 

If  we  assume  that  the  distribution  function  in  gas  is  isotropic,  the  ion  flux  every¬ 
where  in  gas  can  be  derived  as, 


Here,  c  is  the  thermal  velocity  for  a  Maxwell  distribution, 


c  = 


8  kT 


nm-i 


where  rrii  is  mass  of  ion. 

Therefore  the  ion  flux  of  equation  3.14  is, 


At  the  same  time,  from  equation  3.13, 


(3.14) 


(3.15) 


(3.16) 
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Qi,s 

(3.17) 

&  Qa.v 

e  qi,8 

(3.18) 

where  a  is  surface  charge  of  the  liquid  surface  and  e  is  unit  charge. 


The  partition  function  of  the  liquid  surface  is, 
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where  Q\nt  is  the  internal  partition  function  of  ions  in  the  liquid  phase  which 
may  include  rotational,  vibrational  and  excitational  degrees  of  freedom  and  qf™ns  is 
translational  partition  function  of  the  surface  per  unit  area,  h  is  Planck’s  constant. 

The  partition  function  of  gas  can  be  obtained  in  the  same  way,  but  note  that  the 
energy  reference  is  assumed  at  the  liquid  surface, 
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(3.25) 


Qg,v  — 


where  we  define  the  energy  reference  difference  between  liquid  and  gas  as  A  Go 
which  is  called  the  free  enthalpy  or  the  activation  energy,  at  the  same  time,  this  is 
the  energy  barrier  to  cross  over  the  border  between  liquid  and  gas  state.  [38]  Here 
the  system  is  under  equilibrium  with  constant  pressure,  temperature  and  number  of 
molecules  overall  and  this  indicates  that  the  free  energies  do  not  change. 

Therefore  together  with  equation  3.16,  3.18,  3.24  and  3.27,  the  current  density 
becomes, 
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(3.27) 


kT  AGn 
j  =  cr—e  kT 
h 


(3.28) 


This  current  density  can  also  be  considered  as  the  flow  from  the  liquid  surface. 
We  can  apply  this  argument  to  the  liquid  surface  under  electric  field  assuming  that 
the  liquid  surface  is  not  affected  much  by  gas  phase  especially  with  the  field.  [24]  In 
this  case,  the  image  charge  model  3-4  which  is  similar  theory  to  Schottky  model  of 
electron  emission  is  applied. 

Taking  the  x  axis  as  in  figure  3-4  and  applying  an  electric  field  along  it,  the  work 
required  (potential  energy)  to  move  an  extracted  ion  of  charge  e  from  x  to  infinity  is, 
assuming  that  the  binding  force  arises  from  its  image  charge, 


W  = 


(3.29) 


Assuming  the  dielectric  constant  of  material  is  low  enough,  the  force  on  the  charge 
by  its  image  is, 
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Figure  3-4:  Schematic  of  image  charge  model. 
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Here  we  can  see  that  the  image  charge  pulls  back  the  ion  against  the  electric  field. 
The  equation  3.29  becomes, 


W  = 


167re0x 
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(3.31) 


This  energy  is  maximum  when 
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(3.33) 


This  is  called  the  Schottky  depression  which  reduces  the  energy  barrier  by  the 
effect  of  charge  and  field.  Thus,  the  equation  3.28  can  be  modified  as, 


kT 

j  =  a—  exp 
h 


(3.34) 


To  overcome  the  energy  barrier  AG0  solely  by  the  electric  field,  the  Schottky 
depression  needs  to  be  roughly  equal  to  AG0.  Therefore  the  electric  field  necessary 


is: 
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This  indicates  that  the  field  evaporation  occurs  at  around  1-2  V/nm  [26]  which 
corresponds  to  AG'o  of  1.2  -  1.7  eV. 
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Chapter  4 


Ionic  Liquids 

4.1  General  Properties 

Ionic  liquids,  also  called  room  temperature  molten  salts  (RTMS),  are  matter  solely 
composed  by  ions,  anions  and  cations,  and  their  volatility  is  significantly  lower  than 
other  liquids.  [39], [40]  This  property  is  desirable  for  electrospray  thrusters  under  vac¬ 
uum  conditions.  The  volatility  of  ionic  liquids  is  low  because  the  interaction  between 
ions  are  dominated  by  Coulomb- type  forces.  It  is  possible  to  evaporate  ions  from  an 
ionic  liquid  when  an  electric  field  of  enough  strength  is  applied  perpendicular  to  the 
liquid  surface  which  is  shown  in  chapter  3.  The  volatility  limit  of  ionic  liquids  has 
been  recently  probed  spectroscopically  [41]  and  by  distillation  [39].  The  general  con¬ 
clusion  is  that  the  vapor  pressure  of  ionic  liquids,  at  moderate  temperatures,  remains 
negligible.  Another  attractive  point  is  their  high  conductivity  which  enables  a  surface 
reaction  to  strong  electric  fields.  In  addition,  they  are  not  flammable  and  do  rarely 
react  with  other  materials,  thus  avoiding  contamination  of  spacecraft  systems  in  vac¬ 
uum.  However,  some  ionic  liquids  are  sensitive  to  moisture  and  can  be  contaminated 
when  exposed  to  atmosphere  and  also  care  is  needed  due  to  potential  toxicity.  The 
general  characteristics  of  organic  ionic  liquids  are  shown  in  table  4.1.  [42] 

Ionic  liquids  have  found  applicability  in  many  areas  such  as:  green  chemistry  in¬ 
dustry  for  their  stable  and  organic  nature,  gas  handling  using  ionic  liquids  as  solvents 
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Table  4.1:  Typical  Characteristic  of  Organic  Ionic  Liquids 


Low  melting  point 

Treated  as  liquid  at  ambient  temperature 
Wide  usable  temperature  range 

Non- volatility 

Thermal  stability 

N  onflammability 

Composed  by  ions 

High  ion  density 

High  ion  conductivity 

Organic  ions 

Various  kinds  of  salts 

Designable 

Unlimited  combination 

for  wide  variety  of  compounds  and  gases  [43],  nuclear  industry  for  recovery  of  ura¬ 
nium  and  other  metals  from  spent  nuclear  fuel  and  other  sources  [44],  solar  energy 
for  use  as  a  heat  transfer  and  storage  medium  [45]  and  batteries  [46] .  More  recently, 
magnetic  ionic  liquids  have  also  been  reported.  [47] 


4.2  EMI-BF4 

l-ethyl-3-methylimidazolium  tetrafluoroborate  (EMI-BF4)  is  one  of  dialkyl-imidazolium 
cation  family  of  ionic  liquid.  The  general  properties  are  shown  in  table  4.2. 

EMI-BF4  has  one  pair  of  large  anion  and  small  cation  as  shown  in  figure  4-1. 
The  off-balance  of  size  (e.g.  large  cation)  has  been  considered  as  the  reason  why  the 
liquid  state  is  the  most  stable  in  ionic  liquids.  Another  distinguishing  property  of 
EMI-BF4  includes  a  strong  tendency  to  supercool.  [48]  The  electrical  conductivity 
increases  and  the  kinematic  viscosity  decreases  with  increasing  temperature,  same  as 
typical  organic  fluids.  [49]  Another  interesting  characteristic  has  been  identified:  a 


Table  4.2:  EMI-BF4  general  properties 


Melting  point 

[deg  C] 

12.0  -  12.5 

Conductivity  at  room  temp. 

[Si/m] 

1.3 

Dynamic  viscosity  at  room  temp. 

[Pa  .  s] 

0.043 

Density  at  room  temp. 

[kg/m3] 

1517 

Surface  tension  at  room  temp. 

[dyn/cm] 
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phase-change  in  the  vicinity  of  333  K.  [50]  This  is  investigated  by  measuring  diffusion 
coefficients  and  suggests  the  transformation  of  the  diffusion  particle  from  ”  discrete 
ion-pair”  to  ’’individual  ion”  at  temperatures  above  335K  due  to  decomposition  of 
the  EMI-BF4  ion  pair. 

Figure  4-2  shows  partial  charges  of  EMI-BF4  according  to  de  Andrade  et  al.  [51]. 
The  signs  on  atoms  follows  AMBER  atom  type.  AMBER  is  a  type  of  force  field  which 
is  widely  used  for  organic  materials.  The  electric  charge  distribution  of  the  ions  by 
atomic  point  charges  was  calculated  by  quantum  mechanics  (QM)  calculations  on  the 
ab  initio  level. 

As  for  application  to  electrosprays,  from  the  previous  work[52],  it  is  shown  that 
EMI-BF4  is  one  of  the  most  widely  used  liquids  in  electrospray  thruster  research. 
Thrust  efficiency  with  EMI-BF4  is  estimated  to  be  higher  compared  to  other  kinds 
of  ionic  liquids  in  the  work  although  efficiency  is  greatly  affected  by  fragmentation  of 
solvated  ions.  The  detail  of  solvated  ion  is  explained  in  next  section. 
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Figure  4-2:  (a),(b)  Chemical  structural  formulae  with  atom  types  and  partial  charges 
of  EMI+  anion  and  BF4-  cation. 
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4.3  Solvated  Ions 


Solvated  ions  are  defined  here  as  ions  attached  to  neutral  pairs.  For  example,  EMI- 
13  F4  can  have  (EMI-BF4)n  EMI+  for  a  positive  solvated  ion  and  (EMI-BF4)n  BFj  for 
a  negative  solvated  ion.  n  is  defined  as  the  number  of  solvation  which  is  the  number 
of  neutrals  attached  to  the  ion. 

In  general,  the  solvation  is  caused  by  different  types  of  intermolecular  interac¬ 
tions:  hydrogen  bonding,  ion-dipole  and  dipole-dipole  attractions  or  van  der  Waals 
forces.  [53]  Ionic  liquids  are  composed  by  polar  molecules,  thus,  any  of  the  above  may 
happen  to  make  solvated  ions.  In  EMI-BF4  or  in  other  types  of  ionic  liquid,  solvated 
ions  have  been  observed  in  time-of-flight  (TOF)  experiments.  [52], [54], [55], [56], [57] 
Specifically,  it  is  found  by  mass  spectrometry  that  the  species  population  peaks  at 
n  =  0  and  n  =  1,  with  significantly  lower  traces  of  ions  with  n  >  2.  Even  though 
experimental  observations  are  crucial  in  obtaining  information  on  the  macroscopic 
outcome  of  the  field  evaporation  process,  direct  measuring  techniques  are  unable  to 
probe  the  emission  mechanics  at  the  molecular  level.  In  addition  to  this,  it  is  difficult 
to  study  the  stability  of  solvated  ions  through  experiments  alone. 

The  simulation  results  regarding  ion  emission  of  solvated  ions  will  be  presented  in 
chapter  7. 
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Chapter  5 


Numerical  Simulation  Methods 


Continuous  formulations  [58]  an  scaling  analyses  have  been  able  to  describe  in  detail 
the  fluid  dynamics  of  electro  jets  and  the  parametric  relevance  of  ion  emission.  How¬ 
ever,  these  approaches  lose  their  ability  to  track  the  physics  when  the  dimensions  of 
interest  become  similar  in  size  to  the  liquid  molecules.  Molecular  dynamics  (MD)  is 
a  powerful  tool  for  investigating  atomic-scale  phenomena  which  takes  into  account 
intra-  and  inter-  molecular  forces.  This  enables  us  to  segmentalize  most  materials 
into  an  atomistic  scale  at  which  continuum  methods,  like  CFD,  are  not  appropriate. 

In  this  chapter,  the  basics  of  MD  will  be  explained. 


5.1  Molecular  Dynamics  Fundamentals 

Molecular  dynamics  simulations  have  become  a  standard  tool  for  the  investigation  of 
atomistic  problem  such  as  biomolecules  or  structural  failure  [59] .  It  basically  contains 
four  steps  as  shown  in  figure  5-1. 

MD  basically  calculates  the  equation  of  motion  to  obtain  coordinates  and  velocities 
of  next  step.  Forces  in  the  equation  of  motion  are  derived  by  differential  calculus  of 
potentials  (Fj  =  —  dU^i'>)  which  requires  assumptions  depending  on  conditions. 
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Figure  5-1:  Molecular  Dynamics  fundamental  procedure 
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Figure  5-2:  Schematic  of  point  particle  representation  of  covalent  bond 


5.1.1  Potentials  in  MD  Simulations 


MD  does  not  directly  calculate  electrons  and  core  interaction.  They  are  repre¬ 
sented  as  point  particles  and  the  interaction  between  point  particles  includes  repul¬ 
sion  and  attraction  (figure  5-2).  The  black  line  in  right  graph  of  figure  5-2  is  the  sum 
of  energies  of  repulsion  and  attraction  forces  and  particle  dynamics  follow  the  energy 
profile.  The  repulsion  and  attraction  axe  due  to  Pauli’s  exclusion  and  formation  of 
chemical  bond  by  sharing  of  electrons,  respectively. 

Here,  we  assume  the  EMI-BF4  system  has  two  types  of  potentials;  intra-atomic 
potential  and  pair  potential. 
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Energy  U 


Energy  U 


The  harmonic  potential  (figure  5-3)  is  applied  to  the  intra-atomic  potentials  be¬ 
cause  it  can  be  considered  that  atom  bonds  never  break.  As  for  the  pair  potential, 
we  use  electrostatic  (Coulomb)  potential  and  Lenard-Jones  (van  der  Waals)  poten¬ 
tial  [60]. 

The  equation  of  the  potentials  is  expressed  as  follow: 


Kotal  =  E  ~  +  E  K'(d  -  M2 

bonds  angles 

+  E  ^[1  +  cosM-7)]  +  E 

dihedrals  i<j 


Aij  Bij  QiQj 


(5.1) 


where 

=  4eijcr]j 

—  4cijaij 


Here,  is  the  traditional  well-depth  and  is  the  size  parameter.  K  is  the 
energy  coefficient  for  each  kind,  r  is  a  distance  between  atoms  which  are  attached 
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Figure  5-4:  Schematic  of  Lennard- Jones  (LJ)  potential 


each  other  by  covalent  bonds.  9  is  an  angle  within  three  atoms  connected  by  covalent 
bonds,  n  and  7  are  coefficients  for  dihedral  (torsional)  angles  which  is  formed  on 
two  planes.  The  plane  is  defined  by  three  atoms  connected  by  covalent  bonds,  cj) 
is  a  dihedral  angle.  R  is  the  distance  between  pair  atoms.  e0  is  the  permittivity  of 
vacuum,  q  is  a  charge  on  atom.  As  for  suffixes,  they  correspond  to  an  equilibrium 
state.  This  typical  AMBER  force  field  expression  and  the  first  three  terms  represent 
covalent  interactions:  bond  stretching,  bending  and  torsion,  respectively.  The  last 
term  provides  pair  potentials  for  atom-atom  interactions  including  Lenard-Jones  (LJ) 
[60]  and  Coulomb  terms. 


The  Lennard-Jones  potential  is  the  empirical  model  for  van  der  Waals  forces  which 
are  attractive  or  repulsion  forces  between  nonbonded  atoms.  It  is  the  fourth  term  of 
equation  5.1.  The  schematic  of  the  potential  is  shown  in  figure  5-4.  The  well  depth 
e  and  a  corresponds  to  those  in  and  Bjj  in  equation  5.1. 
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Figure  5-5:  Thermodynamic  ensemble  schematic 


5.1.2  Thermodynamic  Ensembles 

Another  important  aspect  of  MD  are  the  thermodynamic  ensembles  which  provide 
a  robust  framework  to  calculate  a  range  of  ’’macroscale”  properties.  MD  provides 
snapshots  of  phenomena  which  are  defined  by  extensive  variables,  which  are  propor¬ 
tional  to  volume  and  mass.  However  the  macroscopic  world  is  presented  by  intensive 
properties  such  as  temperature  and  pressure.  Macroscopic  conditions  translate  to  the 
microscopic  system  as  boundary  conditions.  The  distribution  of  microscopic  states  is 
related  to  the  macroscopic  conditions. 

Effective  method  in  MD  derive  from  so-called  Ergodic  hypothesis,  which  suggests 
that  ensemble  (statistical)  averages  are  equal  to  time  averages  and  all  microstates  are 
sampled  with  appropriate  probability  density  over  long  time  scales. 


1 

n; 


E 


A{i)  —  (A)Ens  — 


(A  /Time 


i...Nt 


(5.2) 


In  a  way,  this  is  a  point  in  common  between  Monte  Carlo  (MC)  simulations  and 
MD  simulations.  MC  and  MD  can  be  compared  in  three  ways. 


1.  How  to  obtain  ensemble  average 

MC:  obtains  ensemble  average  by  random  stepping  scheme. 

MD:  obtains  ensemble  average  by  dynamical  solution  of  equation  (time  history). 
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2.  Assumption 

MC:  Acceptance/rejection  algorithm  leads  to  proper  distribution  of  microscopic 
states. 

MD:  Time  average  equals  ensemble  average  (Ergodic  hypothesis) 

3.  Dynamical  information 

MC:  No  dynamical  information  about  system  behavior  (only  equilibrium  pro¬ 
cesses)  -  no  ’’real”  time 

MD:  Existence  of  dynamical  information  about  system 


5.1.3  Definition  of  Temperature 

We  assume  the  classical  mechanics  (kinetic  theory)  definition  of  temperature: 


(5.3) 


where,  Nf  is  the  total  number  of  degrees  of  freedom  (A/  =  3 A7,  A=  Number  of 
particles)  and  kb  and  T  are  Boltzman  constant  and  temperature  in  Kelvin,  respec¬ 
tively.  Here  average  kinetic  energy  per  degree  of  freedom  is  related  to  temperature 
via  Boltzmann’s  constant. 


5.1.4  Time  Step  Selection 

The  idea  is  to  keep  time  step  as  large  as  possible  due  to  computational  efficiency 
but  not  too  large  to  cause  error.  The  rule  of  thumb  is  that  time  step  must  be  much 
smaller  (  100  . . .  1000  times)  than  the  fastest  time-scale  dynamics  in  the  system,  e.g. 
high  frequency  oscillations  of  light  atoms  (e.g.  C-H  bonds,  O-H  bonds).  Typical 
vibrational  frequencies  is  u>o  —  1013 [sec-1],  thus  it  is  common  to  select  a  time  step  of 
10~15sec  =  1/s.  This  was  also  the  limit  of  time  steps  we  found,  and  stable  results 
could  not  be  provided  for  time  steps  larger  than  1  fs. 
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Chapter  6 


Basic  Numerical  Simulation 
Procedure 


This  chapter  explains  three  steps  carried  out  in  this  research  to  apply  an  external 
electric  field  to  an  EMI-BF4  droplet. 

6.1  Modeling  a  Single  Ionic  Liquid  Molecule 

A  three-dimensional  system  has  been  used  in  this  computational  model  to  ensure  the 
mechanical  integrity  of  the  molecule.  The  optimized  atom  coordinates  (Table  6.1) 
for  a  single  EMI-BF4  molecule  were  obtained  by  the  semiempirical  orbital  program 
MOPAC  6  [61].  We  verify  that  the  cartesian  set  of  coordinates  corresponds  to  the 
structure  of  EMI-BF4  indicated  by  Katsyuba  et  al  [62],  The  structure  includes  loca¬ 
tion  of  anion  relative  to  cation.  We  use  a  force  field  which  was  introduced  for  EMI-BF4 
by  de  Andrade  et  al.  [51].  The  force  field  parameters  are  based  on  the  AMBER  force 
field  [63],  which  contains  a  complete  force  field  for  liquid  state  EMI-BF4  based  on 
the  results  from  both  Quantum  Mechanical  (QM)  and  Molecular  Mechanical  (MM) 
simulations.  Figure  4-2  shows  the  chemical  structural  formulae  which  indicate  atom 
types  and  partial  charges  of  EMI-BF4,  as  they  appear  in  the  results  of  de  Andrade 
et  al. 
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Table  6.1:  Single  EMI-BF4  coordinates.  Unit  in  A 


6.2  Equilibration  of  a  Single  Droplet 


As  the  initial  condition,  64  and  125  EMI-BF4  molecules  are  aligned  in  cubic  with 
potentials  described  in  section5.1.1  (6-1).  Each  molecule  has  the  coordinate  obtained 
in  previous  section  6.1 

Five  different  equilibrium  states  of  125  molecules  are  prepared  at  temperature 
intervals  of  50  K  between  250  K  and  450  K  for  applied  electric  field  simulations. 
Additional  simulations  at  300,  323  and  373  K  are  used  in  non-applied  electric  field 
cases.  In  a  typical  simulation,  atoms  in  the  droplet  are  allowed  to  interact  until 
equilibrium  is  reached  using  a  sequence  of  NVE  and  NVT  ensembles.  NYE  ensemble 
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is  also  called  the  microcanonical  ensemble  in  which  the  number  of  moles,  volume  and 
energy  are  kept  constant.  In  a  similar  way,  NVT  ensemble  is  called  the  canonical 
ensemble  in  which  the  number  of  moles,  volume  and  temperature  axe  kept  constant. 

Excess  energy  stored  in  the  initial  random  coordinate  distribution  is  released  in 
an  NVE  ensemble.  The  procedure  has  several  steps.  In  the  first  step  the  system 
is  initialized  at  10  K.  As  the  simulation  progresses,  the  temperature  increases  until 
converging  to  a  higher  value,  in  some  instances  near  700K.  In  the  second  step,  the 
excess  energy  is  released  by  initializing  the  coordinates  with  the  set  obtained  at  the 
end  of  the  first  step  and  adjusting  the  temperature  back  to  10  K.  This  procedure 
is  repeated  until  the  temperature  converges  to  the  target  temperatures  at  every  50 
K  between  250  K  and  450  K.  The  optimization  of  the  structure  is  then  performed 
with  constant  number  of  moles,  volume,  and  temperature  (NVT,  canonical  ensemble) 
using  a  Nose- Hoover  thermostat  [64]  with  a  temperature  fluctuation  of  100  K  at  the 
target  initial  temperature.  These  simulations  are  carried  out  with  a  cutoff  distance 
of  200  nm  for  both  the  Lennard-Jones  and  the  Coulombic  potentials  thus  accounting 
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for  every  atom  in  the  droplet.  The  equilibrium  states  of  the  EMI-BF4  droplet  of  three 
different  temperatures  are  shown  in  figure  6-2.  The  droplet  has  an  ellipsoidal  shape 
at  least  below  300  K  because  of  the  unsymmetrical  shape  of  the  ions. 

The  condition  of  an  EMI-BF4  droplet  for  125  molecules  after  the  equilibration  is 
temperature:  300K,  droplet  diameter:  4  [nm]  and  EMI-BF4  mass:  lll[u]  for  EMI 
and  87  [u]  for  BF4.  From  here,  the  computed  density  is  1.226  [g/cm3],  which  can 
be  compared  against  the  comparing  the  experimental  value  of  1.24  [g/cm3]  at  295 
[K]  165],  [66]. 
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(c)  400K 

Figure  6-2:  Equilibrium  state  of  EMI-BF4  (a)  300K  (b)  350K  (c)  400K 
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6.3  Droplet  Response  to  Electric  Fields 


Table  6.2:  Simulation  variations. 


E  field  [V/nm 

Temperature  [K] 

jj  of  samples 

Droplet 

External 

1. 0-2.0 

250-450 

5 

Internal 

- 

300,  323,  373 

3 

Droplet  +  Tungsten 

External 

3. 0-6.0 

300 

1 

The  electric  field  is  applied  once  the  equilibrium  state  is  obtained  for  the  different 
initial  coordinates  and  temperatures  (Table  6.2).  The  electric  field  is  applied  in  two 
different  ways:  an  external  constant  and  homogeneous  field  and  through  adding  or 
removing  charges  from  the  droplet.  Here  the  second  way  is  more  realistic  because 
electric  fields  are  far  from  uniform  in  ion  emiting  from  Taylor  cones.  Details  of  the 
droplet  systems  are  explained  in  chapter  7  and  8  and  the  droplet  and  tungsten  system 
is  shown  in  chapter  9.  The  external  constant  electric  field  is  applied  in  the  x  direction. 
In  the  molecular  dynamics  simulation,  the  additional  force  term  (F  =  qE )  is  added 
to  the  fundamental  equation  5.1  to  include  applied  electric  fields. 

Examples  of  ion  emission  for  positive  and  negative  sides  are  shown  in  figure  6-3 
and  6-4.  Due  to  electric  field  in  the  x  direction,  positive  ions  are  emitted  to  the  right 
and  negative  ions  is  to  the  left.  Once  ions  are  emitted,  the  droplet  is  charged  so  it 
starts  moving  along  the  electric  field.  Although  the  NVE  ensemble  is  used  for  the 
simulations,  the  energy  from  the  electric  field  is  added  to  the  system  at  every  time 
step.  Thus,  the  droplet  does  not  cool  down  even  when  in  motion. 
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(a)  0  ps 


(b)  51  ps 


(c)  131  ps  (d)  152  ps 

Figure  6-3:  First  ion  emission  in  positive  side 
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E  field  (x  direction) 


E  field  (x  direction) 


(a)  0  ps 


(b)  46  ps 


(c)  86ps  (d)  106  ps 

Figure  6-4:  First  ion  emission  in  negative  side 
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Chapter  7 


Ion  Evaporation  under  Applied 
Electric  Fields 


7.1  Analysis  of  Electrical  Currents 

Once  the  electric  field  is  applied,  the  droplet  is  polarized  in  the  direction  of  the  field 
(rotating  in  some  cases  due  to  the  induced  dipole  of  the  highly  unsymmetrical  droplet 
(figure  7-1))  and  eventually  experiencing  evaporation  of  non-sol vated  and  solvated 
ions. 

The  emitted  ions  are  observed  at  the  edge  of  droplet  along  the  direction  of  the 
electric  field  and  currents  are  obtained  from  the  ion  dynamics.  The  current  is  esti¬ 
mated  by  counting  the  number  of  ions  for  both  non-solvated  and  solvated  ions  that 
cross  the  liquid  surface  in  a  given  time  interval.  The  total  number  of  positive  and 
negative  emitted  ions  for  each  electric  field  is  approximately  60.  Afterwards,  the 
droplet  collapses  and  retains  nothing  of  the  original  form.  Statistics  of  the  emitted 
ions  are  obtained  until  the  droplet  itself  starts  to  break  up  chaotically  and  observation 
of  ion  emission  is  difficult.  Ion  emission  is  observed  mainly  in  the  electric  fields  for 

1.2  V/nm  through  2.0  V/nm.  The  current  profile  as  a  function  of  applied  electric 
fields  are  examined  in  this  electric  field  range  and  is  shown  in  figure  7-3  [76].  In  the 
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E  field 
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Figure  7-1:  Schematic  of  droplet  rotation 


figure,  the  total  current  includes  both  solvated  and  non-sol vated  ions.  Error  bars 
indicate  standard  deviation  of  five  samples.  The  total  current  profile  is  dominated 
by  non-solvated  ions  in  both  positive  and  negative  cases  and  the  fraction  of  solvated 
ions  is  larger  at  smaller  applied  E  fields.  The  maximum  number  of  solvation  is  n  =  4 
in  the  positive  side  and  n  =  5  in  the  negative  side. 

We  also  examined  current  profiles  as  a  function  of  temperature  to  estimate  the 
activation  energy  of  ion  evaporation  using  Schottky ’s  model  [24] .  The  detail  of  the  ac¬ 
tivation  energy  analysis  is  discussed  in  section  7.4.1.  The  current  profile  is  specifically 
measured  at  an  applied  electric  field  of  1.4  V/nm  ( figure  7-3).  As  expected,  the  cur¬ 
rent  increases  with  temperature  and  the  total  current  is  dominated  by  non-solvated 
ions. 
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Electric  field  [V/nm] 

(a)  (EMI-BF4)„EMI+ 


Electric  field  [V/nm] 

(b)  (emi-bf4)„bf4- 


Figure  7-2:  Current  profile  for  (a)  positive  and  (b)  negative  sides  in  various  applied 
E  fields. 
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Temperature  [K] 

(a)  (EMI-BF4)„EMI+ 


Temperature  [K] 

(b)  (EMI-BF4)„BF4- 


Figure  7-3:  Current  profile  for  (a)  positive  and  (b)  negative  sides  in  various  temper¬ 
ature  with  electric  field  1.4  V/nm. 
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7.2  Local  electric  fields 


Local  electric  fields  are  different  from  applied  electric  fields  due  to  space  charge  dis¬ 
tribution  in  the  droplet.  In  order  to  investigate  the  difference,  the  local  electric  field 
experienced  at  the  position  at  an  extracted  ion  is  estimated  through  the  superposition 
of  the  contribution  of  each  individual  point  charge  and  the  applied  electric  fields. 


and 


E  =  E, 


spacecharge 


+  E, 


external 


(7.1) 


V  •  E  =  —  (72) 

Co 

Example  of  local  electric  fields  along  x,  y  and  z  directions  are  shown  in  figure  7-4 
and  7-5,  respectively.  In  this  case,  the  externally  applied  field  is  directed  along  x 
only.  The  local  electric  field  is  tracked  from  the  beginning  of  the  simulation  until 
the  followed  ion  is  emitted  and  leaves  far  enough  until  the  electric  fields  from  charges 
in  the  droplet  become  negligible.  The  center  of  mass  of  the  target  ion  is  used  for 
the  calculation  and  all  partial  charges  of  every  atom  are  considered,  except  those 
of  the  followed  ion.  Before  emission  occurs,  the  ion  is  exposed  to  strong  negative 
electric  fields  from  its  surroundings  and  is  kept  inside  the  droplet  while  the  field 
magnitude  decreases  throughout  the  relaxation  process.  As  for  ion  emission,  in  a 
macroscopic  sense,  it  is  not  possible  to  measure  such  localized  phenomena.  Ions  leave 
from  a  droplet  when  the  local  electric  field  in  the  x  direction  reaches  its  highest  value, 
while  highest  values  along  y  and  z  do  not  always  happen  when  Ex  is  largest.  The 
Schottky  model  relies  on  a  statistical  ensemble  and  therefore  applies  to  macroscopic 
phenomena.  The  microscopic  details  described  here  are  very  peculiar  to  the  atomistic 
analysis. 
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(a)  x  direction  (b)  y  direction 


(c)  z  direction 

Figure  7-4:  Local  electric  fields  on  a  positive  emitted  ion  in  three  dimension  at  the 
applied  electric  field  1.4  V/nm. 
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Figure  7-5:  Local  electric  fields  on  a  negative  emitted  ion  in  three  dimension  at  the 
applied  electric  field  1.4  V/nm. 
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Figure  7-6:  Local  electric  field  vs.  applied  electric  field  at  the  point  of  ion  emission 
in  300K 


To  investigate  the  trends  of  local  electric  fields  at  the  point  of  ion  emission  in  terms 
of  various  applied  electric  fields,  the  local  electric  fields  are  taken  as  the  average  of 
five  samples  obtained  by  root-mean-square  of  the  electric  field  magnitude  in  all  three 
x ,  y  and  2  directions.  The  local  electric  fields  are  compared  to  the  applied  electric 
fields  as  shown  in  figure  7-6.  Local  electric  fields  are  smaller  than  applied  electric 
fields  for  both  positive  and  negative  ions.  Though  the  number  of  extracted  ions  is 
similar  between  positive  and  negative  sides,  it  can  be  seen  that  positive  ions  (EMI) 
locally  require  higher  electric  fields  probably  because  of  the  unsymmetrical  shape  of 
the  cation.  In  the  experiment  by  Hogan  et  al.  [67],  the  activatoin  energies  differ  by 
approximately  0.25eV. 


62 


Figure  7-7:  Schematic  of  sphere  on  cone  model 


7.3  Estimation  of  propulsive  properties 

Propulsive  properties  of  electrospray  thrusters  are  estimated  using  the  current  and 
electric  fields  relationship  shown  in  figure  7-2  [68]. 

We  assume  the  existence  of  a  Taylor  cone  [1]  with  ion  emissions  issuing  from  its 
apex.  The  electric  field  distribution  along  a  Taylor  cone  surface  is  inversely  pro¬ 
portional  to  \fr,  where  r  is  the  distance  from  the  tip  along  the  surface.  Thus,  ion 
evaporation  could  occur  not  only  from  the  very  tip  of  the  cone,  but  also  fiom  adjacent 
regions  upstream  from  the  tip,  provided  the  field  is  strong  enough.  By  definition,  a 
Taylor  cone  has  infinite  curvature  at  its  apex.  In  real  electrospray  sources,  such  a 
cone  is  an  idealization  which  is  only  accurate  far-  away  from  this  singularity.  Instead 
of  an  ideal  sharp  point,  the  tip  has  some  finite  radius  of  curvature  that  dictates  the 
strength  of  the  field.  The  radius  of  curvature  is  a  function  of  electiostatic  and  liquid 
properties,  such  as  the  applied  voltage,  electrical  conductivity  and  surface  tension. 
To  account  for  the  rounded  apex  at  Taylor  cone,  the  electric  distribution  on  the  liquid 
surface  is  obtained  from  the  sphere  on  cone  (SOC)  model  [69], [55]. 
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In  the  SOC  model,  the  equipotential  lines  are  found  by  solving  the  Laplaces  equa¬ 
tion,  V24>  =  0  in  spherical  coordinates.  The  potential  distribution  is; 


Here,  Pv  is  the  Legendre  function  of  order  u.  Prom  this  potential  distribution,  the 
electric  field  is  derived, 
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Here,  rm  =  fa  is  the  radius  of  the  meniscus,  where  /  is  a  factor  that  determines 
the  surface  equipotential  and  a  correspondes  to  the  model  sphere  radius.  It  can  also 
see  that  the  applied  voltage  is  necessary  to  obtain  the  surface  electric  fields. 

Equation  7.8  is  used  for  the  Legendre  polynomial  considering  the  limit  of  our 
target  x  and  v  (-!<  cosd  <1  and  0<  u  <1). 


PM  =  £(-i)' 


(u  —  n  +  1) . . .  v . . .  (v  +  n)  (  1  —  x 


n= 0 


(n\f 


[v  >  0,  —  1  <  x  <  1) 


(7.11) 

(712) 


Once  the  equipotential  lines  are  obtained  for  a  selected  applied  voltage  between 
the  tip  and  extractor,  the  equipotential  line  whose  angle  relative  to  the  x  axis  is 
approximately  49  degrees  is  selected  as  the  Taylor  cone  surface  model.  Then  the 
electric  field  on  such  equipotential  line  is  found  by  equation  7.5  and  used  for  the 
distribution  along  the  cone  surface.  This  electric  field  distribution  dictates  the  level 
of  ion  evaporation  from  the  Taylor  cone  surface  and  establishes  the  individual  ion 
statistics  for  each  electric  field. 

In  estimating  the  thrust,  velocities  are  measured  when  the  emitted  ions  are  20- 
40  nm  far  from  the  droplet  in  the  direction  of  the  electric  field.  This  distance  is 
not  arbitrarily  chosen.  In  a  true  electrospray  system,  the  field  amplitude  drops  off 
dramatically  with  distance,  yielding  an  equally  pronounced  acceleration  amplitude 
decay.  In  the  MD  simulation  system,  the  ions  are  accelerated  by  constant  electric 
fields.  Here,  we  assume  that  this  20-40  nm  artificial  distance  is  equivalent  to  the 
approximate  terminal  velocities  in  its  physical  correlate,  even  though  weaker  field 
keeps  accelerating  the  ions  through  longer  distances  in  a  real  system.  The  solvated 
ion  break  up  is  taken  into  account  for  the  thrust  on  both  the  positive  and  negative 
sides.  Because  some  of  the  solvated  ions  appear  to  be  in  a  metastable  state,  breakup 
of  solvated  ions  cause  fewer  accelerated  ions.  [54]  Once  the  solvated  ion  breaks  into 
a  neutral(s)  and  a  bare  ion,  only  the  ion  is  accelerated.  Neutral  ion  emissions  are 


65 


Figure  7-8:  Mass  flow  from  droplet 


not  included  in  the  thrust.  The  thrust  density  is  given  as  the  thrust  divided  by  one 
half  of  the  surface  area  of  the  droplet.  The  area  is  approximately  27  nm2.  The  mass 
flow  from  the  droplet  is  measured  and  related  to  total  mass  flow  from  the  surface 
of  the  Taylor  cone.  Data  from  five  different  initial  conditions  are  averaged  out  for 
electric  fields  from  1.2  to  2.0  V/nm.  In  previous  studies,  it  was  found  that  solvated 
ion  emission  is  rarely  observed  above  2.0  V/nm  and  ion  emission  does  not  occur  below 
1.1  V/nm,  at  least  within  the  simulation  time.  [76] 

Based  on  the  ion  evaporation  analysis  using  the  MD  simulation,  the  number  of 
emitted  ions  in  each  electric  field  is  between  20  and  40  in  both  positive  and  negative 
polarities.  After  ions  have  been  emitted  by  this  amount,  the  droplet  itself  loses  its 
consistency  and  breaks  up. 

The  thrust  from  the  droplet  is  shown  in  Figure  7-9.  The  thrust  from  (EMI- 
BF4)nBF4_  is  smaller  than  (EMI-BF4)nEMI+  because  of  the  mass  difference.  As 
expected,  thrust  increases  with  stronger  electric  fields.  Mass  flow  from  the  droplet  is 
shown  in  Figure  7-8.  The  mass  flow  of  positive  ions  is  relatively  random  compared 
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Figure  7-9:  Thrust  from  droplet 


to  the  negative  polarity.  This  tendency  can  be  seen  even  considering  ion  break-up 
because  the  fraction  of  solvated  ions  is  small  when  the  ion  evaporations  occur  at 
more  than  1.3  V/nm.  The  fraction  of  the  solvated  ions  appear  in  Figure  ??.  The 
total  current  is  dominated  by  bare  ions,  and  current  by  solvated  ions  remain  low  even 
under  stronger  electric  fields. 

Figures  7-12,  7-10  and  7-11  show  the  calculated  thrust,  Isp  and  current  from  the 
Taylor  cone  when  2.5  to  3  kV  are  applied  between  the  needle  and  extractor.  These 
values  are  approximately  on  the  same  order  at  same  experimental  results.  [52], [70] 

Because  thrust  is  dominated  by  bare  ions  and  mass  flow  is  not  affected  significantly 
by  the  breakup  of  solvated  ions,  it  can  be  said  the  results  reflect  an  emission  of 
mainly  intact  solvated  ions  even  in  strong  constant  electric  fields.  The  reason  why 
the  breakup  is  hard  to  occur  and  why  the  number  of  emitted  solvated  ions  is  small 
might  relate  to  the  force  field  difference  between  the  stationary  state  and  under  strong 
electric  fields.  A  more  accurate  quantum  model  would  be  required  to  dynamically 
characterize  the  charge  distributions  of  the  emitted  ions.  Such  approach  would  be 
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important  to  establish  the  metastability  of  solvated  ions. 
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Figure  7-12:  Estimated  thrust  from  Taylor  cone 
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7.4  Activation  energy  analysis 


7.4.1  Schott ky  model 

The  activation  energy,  or  the  free  energy  of  the  field  evaporation  process  is  estimated 
using  a  Schottky-type  model  [24],  In  equilibrium,  the  field-enhanced  thermionic  emis¬ 
sion  of  charged  particles  is  described  by, 

-  Mil 

where  a ,  k,  T,  h  and  A  Go  are  the  surface  charge  density,  Boltzmann  constant, 
temperature,  Planck’s  constant  and  activation  energy,  respectively.  The  second  term 
inside  the  exponential  is  commonly  known  as  the  Schottky  depression  and  represents 
the  amount  by  which  the  surface-ion  interaction  barrier  is  decreased  due  to  the  pres¬ 
ence  of  an  electric  field.  This  expression  is  derived  from  statistical  physics  arguments 
in  equilibrium  and  is  therefore  valid  for  aggregate  systems  formed  by  a  large  number 
of  particles.  Nevertheless,  in  the  same  spirit  that  the  statistical  identity  of  a  given 
system  can  begin  to  be  observed  with  a  reduced  number  of  particles  (in  quantum 
statistics,  for  example,  by  finding  the  form  of  the  most  probable  distribution),  the 
nano-droplet  system  composed  of  a  few  thousand  atoms  is  sufficiently  large  for  this 
expression  to  be  applied  with  an  acceptable  degree  of  confidence.  To  obtain  a  numer¬ 
ical  value  of  the  activation  energy,  Equation  7.13  is  compared  against  an  exponential 
fit  of  the  total  emitted  current  density  as  a  function  of  temperature  from  250K  to 
450K  and  a  fixed  electric  field  with  a  magnitude  of  1.4  V/nm  which  is  shown  in 
section  7.1. 

The  exponential  fit  is  used  to  match  the  equation  from  the  Schottky  model  7.13 
to  estimate  the  activation  energy  in  this  section. 

7.4.2  Droplet  energy 

To  take  advantage  of  the  discrete  nature  of  MD  simulations,  the  activation  energy  is 
also  examined  by  keeping  track  of  the  droplet  energy  before  and  after  a  particular  ion 
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is  emitted.  It  is  examined  that  the  temperature  of  droplet  remains  constant  during 
the  emission  process. 

We  assume  the  kinetic  energy  acquired  by  flying  ions  does  not  modify  in  an  appre¬ 
ciable  way  the  internal  energy  of  the  system.  Therefore,  as  an  approximation,  changes 
in  the  droplet  energy  are  only  related  to  changes  in  the  interaction  potentials  within 
the  atoms  that  make  up  the  droplet  and  the  field  evaporated  ions.  The  potential 
energies  include  covalent  bond,  LJ  and  Coulomb  interactions.  The  droplet  energy  is 
expected  to  change  by  a  certain  amount  when  an  ion  is  extracted.  This  energy  change 
is  directly  related  to  the  activation  energy  (or  the  free  energy  for  an  open  system, 
dG\T,p  =  Y2  A4 idNi ,  where  /p  is  the  chemical  potential  and  Nt  the  number  of  particles 
for  species  i),  at  constant  pressure  and  temperature.  This  approach  to  estimate  the 
activation  energy  includes  the  contributions  of  the  different  interaction  potentials  in  a 
direct  way.  Upon  ion  emission,  Coulomb  and  van  der  Waals  forces  practically  vanish 
between  atoms  in  the  droplet  and  the  extracted  ion.  The  extracted  ions  include  both 
non-sol vated  and  solvated  ions. 

Energy  profiles  of  potential  energy  of  droplet  are  obtained  against  time.  The 
approximate  activation  energy  is  derived  from  the  slope  of  energy  vs.  time  character¬ 
istics  and  the  number  of  emitted  ions  per  unit  time.  This  activation  energy  represents 
the  average  energy  change  of  every  ion  emission.  The  analysis  of  activation  energy  is 
done  using  applied  electric  fields  from  1. 2-2.0  V/nm  for  five  different  initial  conditions. 
To  obtain  the  activation  energy,  the  potential  energy  of  the  droplet  tracked  during 
ion  emission.  One  example  of  an  energy  profile  is  shown  in  figure  7-13  This  represents 
the  droplet’s  potential  energy  under  an  electric  field  of  2.0  V/nm  encompassing  the 
extraction  of  four  negative  ion  emission  and  three  positive  ions.  The  activation  energy 
is  approximated  by  the  slope  of  the  energy  during  the  emission  time  interval.  How¬ 
ever,  due  to  the  strong  imposed  electric  field,  the  droplet  deforms  with  an  effective 
reduction  of  the  activation  energy.  Thus,  we  took  only  the  first  5-6  ions  to  obtain  the 
activation  energy  to  avoid  strong  deformations.  The  results  are  shown  in  figure  7-14. 
These  data  are  averaged  over  five  different  initial  conditions.  The  experimental  data 
show  that  nano  droplets  of  EMI-BF4  gives  activation  energies  of  1.83  [eV]  and  1.58 
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Figure  7-13:  Droplet  potential  energy  under  electric  field  of  2.0  V/nm 


[eV]  for  negative  ions  and  positive  ions  respectively.  [67]  From  figure  7-14,  activation 
energies  at  1.2,  1.3  and  1.4  [V/nm]  give  1.5  times  higher  values  than  other  cases. 
However  the  other  values  provide  around  1.8  eV  which  is  a  reasonable  value  when 
compared  with  experiments,  though  further  investigation  is  necessary  to  make  sure 
whether  higher  activation  energies  is  real  or  is  introduced  as  an  artifice  of  the  model. 
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Figure  7-14:  Activation  energy  from  droplet  energy 
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Chapter  8 


Ion  Evaporation  by 


Internal  Fields 


8.1  MD  Approach  to  High  Internal  Electric  Fields 

Based  on  the  Rayleigh  limit  argument,  ion  emission  before  coulomb  fission  was  inves¬ 
tigated  using  MD  simulations.  Coulomb  fission  happens  when  charges  find  a  droplet 
force  the  the  balance  between  electrostatic  and  surface  tension  of  the  liquid  near  the 
limit  of  stability. 

The  electrostatic  pressure  for  a  fully  relaxed  charged  sphere  is, 

Pe  =  ~  (8.1) 

where  a  is  the  surface  charge  density  and  e0  is  the  permittivity  of  vacuum. 

The  surface  tension  pressure  of  the  liquid  droplet  in  vacuum  with  the  spherical 
radius  R  is, 


Pr  =  ^  (8-2) 

where  7  is  the  surface  tension. 

Rayleigh  instability  occurs  when  the  values  of  equation  8.1  and  8.2  become  equiv¬ 
alent. 
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IS, 


qR  =  \J 8'n-2e0'yd 3  (8.3) 

The  electric  field  strength  at  the  droplet  surface  in  the  absence  of  an  external  field 


F_  g 

n  e0d2 

and  the  critical  field  is  from  equation  8.3  and  8.4, 


(8.4) 
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This  indicates  that  the  Coulomb  fission  occurs  under  certain  internal  electric  field 
and  the  electric  field  for  the  fission  has  to  be  larger  than  the  E  field  for  ion  evaporation. 

To  investigate  the  energy  required  for  the  ion  evaporation,  we  model  this  condition 
avoiding  the  fission  using  MD  simulation  removing  ions  from  a  nano-droplet  and 
therefore  forcing  droplet  charging.  The  droplet  contains  64  EMI-BF4  molecules  and 
BF4  ions  are  removed  one  by  one  until  ion  emission  is  observed.  To  achieve  this  we: 


1.  Equilibrate  the  droplet  following  the  same  steps  shown  in  section  6.2  for  three 
different  initial  conditions  at  each  temperature  (300,  323  and  272  [K]). 

2.  Remove  one  BF4  ion  from  the  last  step  of  the  previous  simulation  and  run  a  new 
simulation  using  the  NVT  ensemble  until  the  new  droplet  reaches  equilibrium. 

3.  Check  for  ion  emission.  If  the  emission  is  observed,  stop  the  simulation  and 
restart  the  last  simulation  using  the  NVE  ensemble.  If  ion  emission  is  not 
observed,  repeat  from  step  two. 


Equilibrium  is  checked  by  the  stabilization  of  the  total  energy,  which  is  sum  of 
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energy  of  pair  and  intra  molecular  potentials.  As  for  visualization  of  the  MD  result, 
Visual  Molecular  Dynamics  (VMD)  [71], [72]  is  used  for  the  ion  emission  observation. 


8.2  Activation  Energy  Analysis  as  a  Function  of 
Temperature 

One  case  of  simulation  results  is  shown  in  figure  8-1.  Steps  represent  ions  removed, 
with  the  energy  difference  mainly  caused  by  breakup  of  electrostatic  bonds. 

To  verify  this  point,  we  estimated  the  energy  reduction  per  unit  volume  of  a 
droplet  when  a  single  ion  is  removed.  [37] 

Energy  1  /  e 

Volume  2  0  \47reo-R2 

using  the  volume  of  a  sphere, 

Energv  =  (8'8) 
converting  unit  to  [eV], 

Energy[eV]  =  ^  —  ^  «  0.24[eV]  (8.9) 

€q 7T  il 

Now,  converting  the  energy  unit  in  figure  8-1,  using  100  [g/mol]  as  an  average 
mass  of  EMI  (111  [g/mol])  and  BF4  (87  [g/mol]),  the  energy  step  (approximately  150 
kcal/mol)  is  0.197  [eV],  which  is  close  to  the  value  in  equation  8.9. 


Table  8.1:  Number  of  ions  removed  when  ion  emission  occurs 


Initial  condition  types 

373K 

ver.l 

8 

7 

7 

ver.2 

7 

7 

7 

ver.3 

8 

8 

7 

The  table  9.1  shows  the  number  of  ions  removed  when  ion  emission  occurs.  For 
all  temperatures,  ion  emission  is  observed  when  7  or  8  ions  are  removed.  Although 
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the  energy  profile  also  includes  the  emitted  ion’s  kinetic  and  potential  energy,  the 
large  energy  drop  represents  that  the  pair  potential,  specifically  the  Coulomb  force, 
is  strong  enough  compared  to  the  emitted  ion  energies.  These  emissions  are  not  only 
for  single  ion  emission,  but  also  for  solvated  ion  emission  and  even  droplet  breakup 
as  shown  in  figure  8-2.  Here,  the  energy  drop  for  each  emission  type  takes  almost  the 
same  value  (150  kcal/mol). 

In  order  to  investigate  ion  emission  without  a  thermostat,  NVE  simulations  are 
conducted  for  the  sequence  at  which  ion  emissions  are  observed  in  the  NVT  simula¬ 
tions.  Here,  the  total  energy  is  nearly  constant,  therefore  pair  potential  energies  are 
checked  for  the  energy  profile.  Ion  emission  exist  in  all  initial  conditions/temperatures 
and  we  could  see  energy  drops  for  ion  emission  similar  to  the  NVT  cases  (figure  8-4). 
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Figure  8-1:  Total  energy  profiles  for  ion  removed  simulations  (Initial  condition  version 
1.) 
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(a)  Droplet  breakup Ver.2  373K:  7  BF4  re-  (b)  Single  ion  emission,  Ver.3  300K:  8  BF4  re¬ 
moved  moved 


(c)  Solvated  ion  emission,  Ver.l  300K:  8  BF4 
removed 


Figure  8-2:  Ion  emission  types 
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(a)  Droplet  breakup Ver.2  373K:  7  BF4  re-  (b)  Single  ion  emission,  Ver.3  300K:  8  BF4  re¬ 
moved  moved 


(c)  Solvated  ion  emission,  Ver.l  300K:  8  BF4 
removed 


Figure  8-3:  Total  energy  profile  at  the  sequence  when  ion  emission  is  observed.  Each 
figure  corresponds  to  figure  8-2. 


81 


-3650 


-3900 


100  200  300  400  500 

Time  [ps] 


Figure  8-4:  Pair  potential  energy  profile  of  ion  emission  in  NVE  ensemble:  ver.l,  8 
BF4  removed  at  323K 
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Chapter  9 


Interaction  between  Solid  Plate 
and  Liquid  Droplet 

9.1  Computational  Model 

One  final  model  includes  an  EMI-BF4  droplet  and  a  tungsten  plate  (figure  9-1).  This 
flat  tungsten  plate  represents  the  preliminary  simplified  model  for  the  electrospray 
emitter.  The  size  of  the  tungsten  plate  is  changed  depending  on  the  total  number 
of  molecules.  The  initial  coordinates  for  the  simulations  with  high  electric  fields  are 
obtained  from  the  equilibrium  state  at  near  room  temperature  T  =  300  K.  The  details 
of  the  equilibrium  process  are  discussed  in  section  9.4  and  9.5. 

The  force  field  of  EMI-BF4  is  same  as  before  and  for  the  tungsten  plate,  we 
assume  the  interactions  are  of  the  Lenard-Jones  (LJ)  type,  with  a  Body  Centered 
Cubic  (BCC)  configuration.  For  the  tungsten  atoms,  we  use  the  force  field  parameters 
suggested  by  Tanaka  et  al.  [73] . 


83 


(a)  Atomistic  model  of  tungsten  and  EMI-BF4  droplet 

Nano-droplet 


<J)5  -  20  nm 


(b)  Dimension  of  model  with  the  tungsten  and  8  mol  EMI-BF4 
droplet 

Figure  9-1:  3D  structure  of  the  model  with  tungsten  and  EMI-BF4  droplet 
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9.2  MD  simulation  methodology 

The  following  are  the  steps  for  used  for  the  MD  simulation. 

1.  Equilibration  of  system 

(a)  Equilibration  of  a  nano-droplet  of  EMI-BF4  with  NVT  ensemble. 

(b)  Equilibration  of  tungsten  and  the  nano-droplet  with  NVE  ensemble. 

2.  Simulate  the  system  under  electric  field  conditions. 

•  Apply  constant  and  uniform  electric  field  to  the  system. 

•  Apply  surface  charge  on  the  top  layer  of  the  tungsten  plate. 

For  the  first  step  1(a),  EMI-BF4  molecules  are  placed  arbitrary  in  a  cubic  con¬ 
figuration  as  an  initial  condition.  Once  the  MD  simulation  is  started,  the  cubic 
configuration  gradually  transforms  to  a  spherical  shape.  The  diameter  of  the  droplet 
depends  on  the  number  of  EMI-BF4  molecules.  It  is  approximately  1  nm  to  2  nm 
for  8  and  27  EMI-BF4  molecules,  respectively.  This  simulation  is  conducted  with  the 
NVT  ensemble  for  300  K  using  the  Nose-Hoover  thermostat  [64]  with  a  temperature 
damping  of  100  K.  The  cutoff  distances  are  8  A  and  100  A  for  LJ  and  Coulombic  po¬ 
tentials,  respectively.  The  time  step  is  1  fs  and  the  maximum  total  simulation  time  is 
50  ns  which  corresponds  to  50  million  iterations.  The  periodic  boundary  condition  is 
not  used  because  this  is  an  isolated  floating  system.  The  completion  of  equilibration 
is  checked  both  by  visualization  and  energy  profiles.  After  the  equilibrium  is  reached 
with  the  NVT  simulation,  the  final  trajectory  of  1(a)  is  used  for  the  initial  condition 
of  1(b).  Here,  the  tungsten  plate  is  added  to  the  system  and  the  distance  between 
the  nano-droplet  and  the  plate  is  set  to  approximately  3  A.  The  initial  condition  of 
the  tungsten  plate  has  BCC  configuration  with  a  lattice  constant  of  3.165  A.  This 
tungsten  plate  initial  configuration  is  obtained  by  several  trials  using  different  shapes. 
The  details  are  in  section  9.5. 

Once  the  equilibration  is  achieved  in  step  1(b),  the  tungsten  and  the  nano-droplet 
are  exposed  to  electric  fields.  We  investigated  two  types  of  electric  field  circumstances 
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as  mentioned  above.  The  former  model  provides  a  simpler  configuration  and  has  been 
discussed  in  several  studies.  [74]’  [75]’  [76]  However,  because  it  is  only  consistent  with 
a  realistic  model  very  near  to  the  surface,  we  consider  also  a  model  using  surface 
charges,  which  gives  a  decaying  electric  field  which  is  more  realistic.  The  details  of 
the  surface  charge  approach  are  discussed  in  section  9.2.  The  constant  electric  field 
is  applied  perpendicular  to  the  surface  of  the  tungsten  plate.  In  these  calculations, 
the  position  of  all  tungsten  atoms  are  fixed  due  to  the  severe  computational  cost  that 
otherwise  would  arise.  Also,  to  make  the  simulations  faster,  a  neighbor  list  is  applied. 
The  skin  distance  for  the  neighbor  list  is  set  to  10  A.  A  shrink  wrapping  algorithm 
is  used  for  ion  emission  which  have  a  possibility  of  leaving  the  simulation  box.  The 
boundary  conditions  are  not  periodic  in  this  calculation,  either,  and  the  ensemble  is 
NVE  with  a  time  step  of  1  fs.  The  initial  temperature  is  300  K. 


9.3  Surface  Charge  on  Tungsten  Plate 

To  simulate  a  charge  distribution,  charges  are  applied  to  atoms  on  the  top  layer  of 
the  tungsten  plate,  such  that 

=  ^  =  eo  E  (9.1) 

where  asur ,  Q.  S,  e0  and  E  are  surface  charge,  total  charge  on  the  surface,  surface 
area,  permittivity  and  electric  field,  respectively.  The  area  is  obtained  from  the  lattice 
constant  as  shown  in  figure  9-2.  Here, 


Q  =  Nq  (9.2) 

where  N  and  q  are  the  number  of  atoms  at  the  top  layer  of  the  tungsten  plate 
and  charge  on  each  atom,  respectively.  Therefore,  the  charge  for  each  atom  is 


9  =  CE-E  (9.3) 

This  amount  of  charge  is  applied  to  provide  an  electric  field  E  on  the  surface.  The 
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First  layer  of  tungsten  plate 


Tungsten 

plate 


Nano-droplet 
of  EMI-BF4 


a:  lattice  constant 
S:  area  of  one  unit  cell 


Infinite  plate 


Figure  9-2:  Schematic  of  the  first  layer  of  the  tungsten  plate. 

electric  field  is  varied  to  observe  ion  emission  of  EMI-BF4. 

It  is  important  to  point  out  that  a  uniform  charge  distribution  on  the  finite  sim¬ 
ulated  plate  does  not  produce  an  equipotential  surface  and  as  a  consequence  the 
corresponding  field  is  not  necessarily  normal  to  the  surface.  This  is  a  good  approxi¬ 
mation  only  if  the  size  of  the  plate  is  much  larger  than  the  size  of  the  droplet,  which 
is  not  the  case  of  our  simulations,  mostly  because  of  computational  cost.  In  future 
work  we  plan  to  address  this  problem  by  incorporating  periodic  boundary  conditions 
(see  figure  9-2,  left  drawing)  such  that  an  infinite  plate  is  simulated  to  produce  a 
normal  field  from  a  uniform  charge  distribution. 


9.4  Equilibration  of  Tungsten  Plate 

Because  it  is  necessary  to  use  LJ  potential  for  the  tungsten  plate  to  calculate  the 
interactions  under  non-periodic  boundary  conditions,  there  is  a  computational  limi- 
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tation  to  the  thickness  of  the  tungsten  plate.  We  tried  several  patterns  to  reduce  the 
computational  cost.  The  cutoff  distance  is  20  A  for  LJ  potential  and  the  Coulombic 
force  is  not  calculated  because  tungsten  atoms  are  initially  uncharged. 

Figure  9-3  shows  three  examples  of  tungsten  plates.  The  left  figures  (figure  9- 
3  (a),  (d)  and  (g))  show  model  schematics,  and  the  center  (figure  9-3  (b),  (e)  and 
(h))  and  right  (figure  9-3  (c),  (f)  and  (i))  figures  correspond  to  snapshots  of  initial 
coordinates  and  time  varying  coordinates,  respectively.  The  round  model  1  (figure  9-3 
(d),  (e)  and  (f))  is  built  with  a  concept  that  provides  a  thick  tungsten  plate  where 
the  nano-droplet  is  placed  while  other  parts  have  smaller  thickness  to  reduce  the 
number  of  tungsten  atoms.  The  square  model  (figure  9-3  (a),  (b)  and  (c))  and  round 
model  2  (figure  9-3  (g),  (h)  and  (i))  have  uniform  thickness  of  3  lattices.  It  can  be 
seen  that  the  there  are  deformations  in  the  square  model  and  the  round  model  1. 
The  deformation  of  the  square  model  is  not  critical  compared  to  the  round  model  1, 
however  the  number  of  atoms  of  the  round  model  2  are  less  than  the  square  model. 
Also,  the  round  model  2  perfectly  keep  its  configuration  as  shown  in  figure  9-3  (h) 
and  (i).  Therefore  we  use  the  round  model  2  for  further  simulations. 


9.5  Equilibration  of  Tungsten  Plate  and  EMI-BF4 
Nano-droplet 

Equilibration  of  the  tungsten  plate  and  the  EMI-BF4  nano-droplet  is  done  for  two 
different  system  sizes  with  8  and  27  molecules  EMI-BF4  droplet.  The  diameters  of  the 
tungsten  plates  are  8  nm  and  10  nm,  respectively.  The  total  number  of  atoms  are  2499 
and  5343  for  these  systems.  The  combination  of  (a),  (b)  of  figure  9-4  and  9-5  shows 
the  initial  condition  and  the  equilibrated  condition  of  8  and  27  EMI-BF4  molecules, 
respectively.  Compared  to  the  sole  equilibration  of  droplets,  these  equilibrations  need 
a  longer  cutoff  distance  for  LJ  potential  because  of  the  relatively  large  LJ  coefficient 
of  tungsten  atoms.  Here,  we  use  30  A  and  80  A  for  8  and  27  EMI-BF4  molecules 
respectively.  From  figure  9-4  and  9-5,  it  can  be  seen  that  the  equilibration  condition 
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15  nm 


(a)  Schematic  of  square  model.  (b)  Square  model,  0  ps.  (c)  Square  model,  820  ps. 


1  h-»15nm  H  l 

2  lattices  t  '  j  1  lattice 

<t>  5  nm 

(d)  Schematic  of  round  model 

1. 


3  lattices 


(g)  Schematic  of  round  model 

2. 


(e)  Round  model  1,  Ops. 


(h)  Round  model  2,  0  ps. 


(f)  Round  model  1,  1090  ps. 


Figure  9-3:  (a)  (b)  (c)  Square  model  of  tungsten  plate,  (d)  (e)  (f)  Round  model  1  of 
tungsten  plate,  (g)  (h)  (i)  Round  model  2  of  tungsten  plate 
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is  totally  flat  and  both  EMI+  and  BFJ  ions  completely  wetted  the  tungsten  surface 
due  to  a  strong  attraction  by  tungsten  atoms.  For  a  more  precise  investigation  of  ion 
emission  from  a  Taylor  cone,  a  thicker  liquid  layer  might  be  needed. 

Figure  9-6  is  the  total  energy  profile  throughout  the  equilibration  runs.  The  8- 
molecule  system  is  equilibrated  approximately  at  50  ns,  and  the  27-molecule  system 
converge  at  around  30  ns.  The  physical  simulation  time  is  4  days  for  350  ps  in  the  27 
molecules  case,  even  using  a  neighbor  list  and  fixed  tungsten  atoms.  Here,  the  force 
on  tungsten  atoms  are  not  calculated  and  are  set  to  zero.  The  most  expensive  part 
is  the  calculation  of  the  LJ  and  Coulombic  potentials  with  long  cutoff  distances.  To 
study  a  larger  system,  it  will  be  necessary  to  consider  certain  techniques  such  as  the 
multiple  time  step  method[77]. 

9.6  Tungsten  plate  and  nano-droplet  under  elec¬ 
tric  fields 

Here  we  investigate  tungsten  plate  and  nano-droplet  under  two  sources  of  electric 
fields.  Table  9.1  indicates  the  ion  emission  existence  in  four  different  conditions: 

a.  Ion  emission  under  external  constant  electric  field  with  8  EMI-BF4  molecules. 

b.  Ion  emission  with  surface  charge  with  8  EMI-BF4  molecules. 

c.  Ion  emission  under  external  constant  electric  field  with  27  EMI-BF4  molecules. 

d.  Ion  emission  with  surface  charge  with  27  EMI-BF4  molecules. 

Figure  9-7  shows  an  EMI+  ion  emission  from  the  27-molecule  system.  The  electric 
field  is  applied  in  the  x  direction  which  is  perpendicular  to  the  surface.  Ion  emission 
from  an  ionic  liquid  surface  usually  occurs  in  electric  fields  of  1.0  V/nm  -  2.0  V/nm. 
However,  here,  the  electric  field  needed  for  ion  emission  is  higher  than  those  values 
because  ions  are  directly  attached  to  the  metal  surface. 

At  electric  fields  stronger  than  3.5  V/nm,  ion  emission  along  the  x  direction  is 
observed  under  the  external  constant  electric  field  condition,  but  no  axial  ion  emission 
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(a)  Initial  condition  with  8  EMI-BF4  molecules. 


(b)  Equilibrated  condition  with  8  EMI-BF4  molecules. 


Figure  9-4:  Equilibration  of  8  EMI-BF4  molecules  and  tungsten  plate. 
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(a)  Initial  condition  with  27  EMI-BF4  molecules. 


(b)  Equilibrated  condition  with  27  EMI-BF4  molecules. 


Figure  9-5:  Equilibration  of  27  EMI-BF4  molecules  and  tungsten  plate. 
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(a)  8  EMI-BF4  molecules  and  tungsten  plate. 


(b)  27  EMI-BF4  molecules  and  tungsten  plate. 


Figure  9-6:  (a)  (b)  Total  energy  profile  of  equilibration  with  EMI-BF4  molecules  and 
tungsten  plate. 
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Figure  9-7:  EMI+  Ion  emission  from  tungsten  surface. 


is  observed  through  a  surface  charge  distribution.  However,  in  the  system  with  a 
surface  charge,  it  can  be  seen  that  a  single  molecule  is  emitted  from  the  the  edge  of 
the  tungsten  surface  where  the  electric  field  is  strongest.  This  emission  occurs  at  a 
surface  charge  distribution  which  corresponds  to  6.0  V/nm.  This  reinforces  the  notion 
that  a  modified  surface  charge  model  is  needed  to  model  perpendicular  ion  emission 
The  number  of  ions  emitted  is  larger  in  the  system  with  27  molecules  as  shown 
in  table  9.1.  From  the  ion  emission  tendency,  the  LJ  potential  of  tungsten  atoms 
seems  to  work  properly,  though  further  investigation  with  more  layers  of  EMI-BF4  is 
needed  to  compare  against  experimental  values.  Also,  there  is  the  more  interesting 
aspect  of  this  type  of  simulations,  namely,  the  ionic  (electrochemical)  interactions  of 
counterions  with  metal  structure.  That  might  require  some  form  of  reactive  MD. 
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Table  9.1:  Ion  emission  profile  for  27  EMI-BF4  molecules  and  tungsten  plate,  (a) 
Ion  emission  under  external  E  field  with  8  EMI-BF4  molecules,  (b)  Ion  emission  with 
surface  charge  with  8  EMI-BF4  molecules,  (c)  Ion  emission  under  external  E  field 
with  27  EMI-BF4  molecules,  (d)  Ion  emission  with  surface  charge  with  27  EMI-BF4 
molecules.  Parenthesis  indicates  the  time  when  the  first  ion  emission  is  observed. 


Electric  field  [V/nm 

(a) 

(b) 

(c) 

(d) 

3.0 

- 

- 

- 

- 

3.5 

O  (138  ps) 

- 

0  (139,  140  ps) 

- 

4.0 

0  (41,  182  ps) 

- 

O  (19,  92,  113  ps) 

- 

5.0 

0(1,2,  7  ps) 

- 

O  (3,  4,  16,  17  ps) 

- 

6.0 

O  (Every  1  ps) 

O  (198  ps) 

O  (Every  1  ps) 

O  (172  ps) 
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Chapter  10 


Conclusions  and  Future  Work 

10.1  Summary  of  Results  and  Contributions 

Atomistic  modeling  has  the  potential  to  be  a  versatile  tool  to  reveal  the  physics 
of  electrosprays.  Preliminary  studies  for  electrospray  thrusters  have  been  done  by 
molecular  dynamics  simulations  of  ionic  liquid  droplets.  The  goal  is  to  establish  the 
first  milestone  of  multi-scale  modeling  of  electrosprays  to  reveal  their  mechanisms 
mainly  based  on  the  analysis  of  ion  emission  current  and  activation  energy. 

Current  was  measured  by  keeping  track  of  the  number  of  emitted  ions  from  ionic 
liquid  droplets  under  two  conditions:  fixed  temperature  (300  [K])  at  various  electric 
fields  and  fixed  electric  field  (1.4  [V./nm])  at  various  temperatures  (250  -  450  [K]). 
Electric  current  analysis  was  conducted  using  simulation  results  for  the  first  time  and 
the  obtained  total  current  density  was  reasonably  well  matched  with  experimental 
values.  These  results  show  ion  emission  of  both  solvated  ions  and  non-sol vated  ions 
with  a  largest  number  of  solvation  of  n  —  4  in  the  positive  side  and  n  =  5  in  the 
negative  side.  But  n  —  0  was  by  far  more  common.  Activation  energy  analysis  was 
made  by  fitting  current  vs.  temperature  curves  with  a  Schottky-type  model.  This 
resulted  in  slightly  lower  values  for  the  activation  energy  compared  to  experimental 
data  with  an  average  error  of  approximately  10  percent.  The  potential  energy  of 
droplets  is  also  measured  for  various  electric  fields  (1. 2-2.0  [V/nm])  as  an  alternative 
to  study  activation  energy  under  external  electric  fields.  It  is  found  that  this  method 
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is  greatly  affected  by  droplet  deformation  so  only  the  first  five  or  six  emitted  ion  are 
used  for  the  analysis.  This  analysis  agrees  relatively  well  with  the  latest  experimental 
results  for  nano-droplets.  Activation  analysis  through  internal  electric  fields  was  also 
considered.  The  energy  drop  corresponded  well  to  a  rough  analytical  estimation.  The 
summary  of  activation  energy  for  each  method  is  shown  in  table  10.1. 


Table  10.1:  Approach  anc 

activation  energies. 

Efield 

Method 

Activation  energy  per  ion  [eV] 

External 

Schottky  fitting 

1.4 

Droplet  potential  energy 

1.8 

Internal 

NVT 

0.2 

NVE 

0.1 

Propulsion  properties  have  been  investigated  applying  the  ion  emission  results 
from  droplets  and  electric  field  distributions  along  the  liquid  surface  using  the  sphere- 
on-cone  model.  Thrust,  Isp  and  current  from  one  Taylor  cone  were  calculated  for  both 
positive  and  negative  polarities. 

A  final  model  was  presented  for  the  interaction  between  a  liquid  droplet  and  a 
solid  have  been  observed  using  full  AMBER  force  field  model.  Electric  fields  from 
a  surface  charge  distribution  was  applied  which  gives  a  decaying  field  model  that  is 
closer  to  reality.  The  results  indicate  that  the  surface  charge  would  work  but  the 
modification  for  the  charge  distribution  and  force  field  are  needed  for  more  realistic 
model. 

Molecular  dynamics  simulations  yield  reasonable  predictions  for  the  activation 
energy  of  ionic  liquids.  Further  research  is  necessary  as  mentioned  in  section  10.2, 
however,  atomistic  modeling  is  a  promising  way  to  understand  this  nano-scale  phe¬ 
nomena. 
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10.2  Recommendations  for  Future  Work 


10.2.1  MD  Simulation  of  Droplets 

As  indicated  in  chapter  6  and  7,  there  are  several  problems  in  the  droplet  simulations. 

Among  them: 

a.  Moving  droplet 

Under  external  electric  fields,  the  floating  droplet  moves  along  the  electric  field 
when  it  is  charged  due  to  ion  emission.  It  is  possible  to  analyze  ion  emission 
with  the  moving  droplet  but  it  prevents  from  a  simple  kinetic  and  potential 
energy  analysis  and  may  affect  the  energy  calculations  especially  when  a  larger 
droplet  is  used  with  more  ion  emission  because  it  also  gives  shorter  time  to 
accelerate  the  droplet  itself. 

The  suggested  solution  is  fixing  four  or  five  molecules  placed  at  the  center  of 
the  droplet.  This  is  artificial  but  it  should  not  affect  too  much  the  behavior  of 
surface  ions  when  a  droplet  is  large  enough  and  the  atoms  keep  their  partial 
charges  and  LJ  potentials.  The  fixing  should  be  made  after  the  droplet  is  fully 
relaxed  so  that  no  more  rotation  occurs  due  to  dipoles. 

b.  Number  of  molecules 

The  largest  number  of  atoms  in  this  work  is  3000,  but  in  general  MD  is  able  to 
handle  more  than  ten  to  hundred  billions  of  atoms.  However,  this  is  possible 
when  special  treatments,  such  as  periodic  boundary  conditions  or  a  reasonable 
neighbor  list  are  applied  and  the  code  is  run  on  “Teraflop”  computers  which 
require  more  than  1000  Opteron  CPUs  and  large  memory.  (The  estimation 
of  computational  time  can  be  made  with  flops  per  clock  per  CPU  for  each 
computer.)  One  of  the  reasons  for  the  limited  number  of  molecules  is  the  time 
step  restriction  (section  5.1.4)  which  requires  approximately  1,000,000  iteration 
steps  to  cover  only  1  [ns]  trajectories.  In  addition  to  massive  iterations,  pair 
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potentials,  especially  long-range  Coulomb  potentials,  increase  computational 
burden  due  to  the  summations  of  every  combination  of  particles  thus  scaling  as 
N2.  Cutoff  distances  and  neighbor  lists  can  reduce  this  amount  but  this  trade¬ 
offs  with  the  accuracy  of  the  potential  calculation.  Also,  periodic  boundary 
conditions  are  difficult  to  apply  to  droplet  systems  which  axe  isolated  from 
their  surroundings.  Thus,  to  tackle  on  larger  systems,  fast  CPUs  and  large 
memories  or  appropriate  modeling  of  systems  are  necessary.  To  observe  an  ion 
emission  from  liquid  surfaces,  it  might  be  necessary  to  analyze  a  large  plate 
of  liquid  instead  of  droplets,  including  periodic  boundary  conditions  and  apply 
external  electric  fields.  If  it  is  necessary  to  keep  investigating  droplet  physics, 
NSF  supercomputer  is  available  for  larger  computations.  [78] 

c.  Fraction  of  solvated  ions 

As  shown  in  section  7.1,  the  fraction  of  solvated  ions  against  non-solvated  ions 
are  small  although  experimental  results  show  that  monomer  and  dimer  almost 
give  same  percent  of  total  current  in  negative  side  and  4:6  in  positive  side.  [52] 
This  indicates  that  solvated  ions  in  the  simulation  are  easy  to  break,  id  est, 
there  is  a  possibility  that  the  coefficients  of  the  force  field  are  not  appropriate 
in  this  case.  The  force  field  is  validated  in  equilibrium  states  [51],  however, 
apparently  not  under  strong  electric  fields.  In  reality,  the  external  electric  fields 
might  change  the  partial  charge  which  atoms  have  and  even  the  van  der  Walls 
potential.  Another  possibility  is  the  effect  of  different  electric  field  distributions: 
Taylor  cone  and  constant  electric  fields,  a  condition  that  can  be  represented  as 
a  droplet  placed  between  parallel  plate  electrodes.  Further  consideration  is  also 
required  for  the  value  of  dielectric  constant  of  the  ionic  liquid. 

10.2.2  MD  Simulation  of  Liquid  and  Solid 

Further  investigation  is  necessary  mainly  for  the  selection  of  an  adequate  force  field. 

Currently  it  is  impossible  to  mix  up  two  different  force  fields  in  one  MD  simulation 

and  the  AMBER  force  field  is  the  best  option  for  organic  materials  but  not  for  metals. 
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The  reason  that  metals  have  different  force  fields  can  be  described  by  their  bonding 
system.  Organic  materials  have  pair  potentials  as  inter  molecular  potentials,  but  on 
the  other  hand,  metals  have  multi-body  potentials  in  which  bond  strength  is  weaker 
at  the  surface  of  crystals.  This  is  affected  by  quantum  mechanical  effects  that  describe 
the  influence  of  the  electron  gas.  Embedded  atom  method  (EAM)  potentials  are  used 
as  force  fields  for  metals,  however,  it  does  not  have  options  to  calculate  liquids.  It 
is  known  that  the  first-principle  based  ReaxFF  has  an  ability  to  calculate  two-phase 
systems  although  it  is  not  capable  get  to  calculate  ionic  liquids  because  lacks  many 
atomic  interaction  types.  We  might  be  able  to  wait  for  ReaxFF  to  add  new  atoms 
but  there  is  another  option  to  create  a  force  field  called  CVFF  using  appropriate 
software.  [79]  Technically  it  is  possible  to  obtain  the  force  field  by  direct  calculation 
but  this  is  not  recommended  because  the  process  takes  too  much  time  in  order  to 
adjust  the  potential  curve  until  obtaining  the  appropriate  values  of  parameters  for 
every  intra-  and  inter-  atomic  potential.  The  software  for  CVFF  force  field  is  material 
“Material  Studio”  and  the  CVFF  force  field  also  works  on  LAMMPS.  [80] 

Electrochemical  effect  between  solid  and  liquids  would  require  a  quantum  method 
because  molecular  dynamics  simulation  does  not  directly  calculate  the  motion  of 
electrons.  The  problem  of  quantum  methods  is  that  they  handle  smaller  number  of 
atoms  than  MD.  Thus  it  might  be  difficult  to  see  mesoscopic  phenomena  such  as 
Taylor  cone  with  QM  calculation  due  to  the  limited  capacity  and  processing  speed  of 
machines. 

10.2.3  Other  Simulation  Techniques 

Other  than  molecular  dynamics  simulation,  PIC  (Particle  in  Cell)  model  is  one  of 
options  to  investigate  behavior  of  ion  jets.  The  detail  is  described  in  Appendix  C. 
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Appendix  A 


Validation  of  MD  simulation 

A.l  Relaxation  Time 

The  relaxation  time  of  the  EMI-BF4  droplets  under  electric  fields  are  checked.  The 
droplets  have  125  EMI-BF4  molecules  and  these  were  all  equilibrated  at  300  [K]  using 
thermostat  [64]  before  the  electric  fields  are  applied.  The  time  when  first  ion  is  emitted 
is  assumed  when  the  relaxation  is  completed  and  the  time  of  five  different  samples 
in  each  electric  field  strength  are  averaged  out.  Figure  A-l  shows  the  obtained  result 
both  for  positive  and  negative  emissions.  We  can  see  that  the  simulation  values  are 
in  same  order  with  typical  relaxation  time  which  is  in  order  of  100-200  ps  obtained 
from  arguments  in  section  3.1. 
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Figure  A-l:  Charge  relaxation  time 
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A. 2  Evaporation  Temperature 


The  evaporation  temperature  is  investigated  for  EMI-BF4  droplet  with  64  molecules. 

Starting  from  room  temperature,  temperature  is  increasing  using  Nose  Hoover 
thermostat  with  dumping  of  100  [K]  as  shown  in  figure  A-2.  Force  fields  of  molecules 
are  same  as  before  [51],  but  there  is  no  electric  field  in  this  case.  The  neutral  ion 
emission  was  observed  at  approximately  1140K  which  is  higher  than  experimental 
value  (approximately  400  [K]).  However,  the  experimental  value  is  not  obtained  by 
a  nano  droplet  neither  emission  of  one  molecule.  Therefore  further  investigation  is 
needed  to  explore  the  evaporation  temperature  of  ionic  liquids. 


Figure  A-2:  Temperature  transition 
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Appendix  B 


Computation 

B.l  Molecular  Dynamics  Software 

The  LAMMPS  (http://lammps.sandia.gov/)  MD  software  has  been  used  throughout 
this  research.  All  the  descriptions  and  explanations  of  commands  for  input  files  can 
find  on  the  website.  To  make  LAMMPS  run  on  the  SPL  cluster  machine,  OpenMPI 
option  is  selected. 


B.2  Visualization  Software 

B.2.1  General  Information 

Visual  Molecular  Dynamics  (VMD)  was  used  for  visualization  of  MD  simulation  re¬ 
sults.  VMD  is  developed  by  theoretical  and  computational  biophysics  group  in  Uni¬ 
versity  of  Illinois  at  Urbana-Champaign  and  is  available  for  free  from  the  link  below. 

http://www.ks.uiuc.edu/Research/vmd/ 

Instructions  and  tutorials  can  be  also  found  on  the  website.  As  for  visualization 
of  LAMMPS  results,  we  use  dump  and  dcd  files  to  load  on  the  VMD. 


107 


B.2.2  Ion  Tracking 


An  ion  tracking  script  was  created  to  make  movies  following  moving  droplets.  This 
script  was  based  on  Dr.  Axel  Kohmeyer’s  work  and  the  details  can  be  found  in  section 
6.2  of  the  link  below. 

http:  /  /  sites.google.com/site/akohlmey/ redirect/cpmd-vmd.pdf 

The  following  is  the  script  for  our  EMI-BF4  droplet.  Underlines  show  parts  where  we 
needed  to  make  modifications  depending  on  the  files  and  target  ions. 


tt - 

jj  ion  tracking  code  for  EMI-BF4 

ti - 


jj  load  trajectory: 

mol  new  {(dump  file  name)}  type  lammpstrj  waitfor  all 
mol  addfile  {(dcd  file  name)}  type  dcd  waitfor  all 
jjmol  delrep  0  top 

jjmol  selection  type  123456789  10 

set  molid  1 
set  hoffs  0 

j)  let  all  selections  be  recalculated  for  each  frame 

jj  and  smooth  the  trajectory  a  little  bit  for  all  representations 

jj  that’s  part  two  of  the  magic. 

set  n  [molinfo  Smolid  get  numreps] 

for  {set  i  0}  {$i  <  $n}  {incr  i}  { 

mol  selupdate  $i  Smolid  on 
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mol  smoothrep  $molid  $i  2 

} 


proc  do_realign  {args}  { 

global  molid  chist  hcount  hoffs  dhist 

jj  this  is  the  axis  to  align  to 
jjset  avec  [vecnorm  {1.0  0.0  0.0}] 
jj  this  is  the  sliding  window  size 
set  asize  5 

(j  initialize  the  cache  counters 
if  ([info  exists  hcount])  {  }  else  { 
set  hcount  0 
set  hoffs  -1 
} 


jj  find  center 

set  sel  [atomselect  Smolid  “underlineindex  (atom  ID)  or  index  (atom  ID)”] 

lassign  [$sel  get  {x  y  z}]  emil  emi2 

set  cent  [vecscale  [vecadd  erralemi2]  0.5] 

jjset  dir  [vecsub  ermlemi2] 

jj  store  data  in  cache  for  sliding  window  averaring 
if  {$hcount  <  $asize}  then  {  incr  hcount  } 
incr  hoffs 

if  {Shoffs  >=  $hcount}  then  {  set  hoffs  0  } 
set  chist($hoffs)  Scent 
jjset  dhist(Shoffs)  $dir 
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tt  calculate  averages 
set  csum  [veczero] 
j}set  dsum  [veczero] 

for  {set  i  0}  {$i  <  $hcount}  {incr  i}  { 
set  csum  [vecadd  $csum  $chist($i)] 

{[set  dsum  [vecadd  $dsum  $dhist($i)] 

} 

set  csum  [vecscale  [expr  1.0/[expr  Shcount  *  1.0]]  $csum] 
jjset  dsum  [vecnorm  $dsum] 

j]  get  rotation  axis 

{[set  rvec  [vecnorm  [veccross  $dsum  $avec]] 

(J  set  origin 
molinfo  Smolid  set 
{center_matrix} 

[list  [trans  origin  $csum]  ] 

0[trans  axis  $rvec  [expr  acos([vecdot  $dsum  $avec])]  rad]] 

j[  clean  up  selections 
$sel  delete 
} 


scale  to  0.01 

trace  variable  vmd-frame(l)  w  do_realign 


(J  go  back  to  the  start  of  the  trajectory, 
{[animate  style  rock 
[[animate  goto  start 
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B.3  SPL  Cluster  Machine  and  Execution  of  LAMMPS 


All  the  simulation  of  this  work  has  been  done  using  SPL  cluster  machine  (spl.mit.edu). 
It  has  16  nodes  with  four  AMD  Opteron  cores  per  node  and  two  nodes  with  8  Shang¬ 
hai  cores  per  node.  OS  is  Red  hat  enterprise  Linux  and  all  the  basic  options  for 
numerical  simulations  are  available.  However,  it  is  necessary  to  modify  bash  file  of 
an  account  to  install  the  LAMMPS  on  the  account  before  making  LAMMPS.  The 
“.bashrc”  file  at  a  home  directory  of  an  account  needs  to  be: 


D  .bashrc 

jj  Source  global  definitions 
if  [  -f  /etc/bashrc  ];  then 
.  /etc/bashrc 

fi 

jj  User  specific  aliases  and  functions 
PA  T H=$PA  TH: /home /(user  name ) /bin 


Modification  of  a  make  file  is  also  required  depending  on  the  bash  settings.  The 
“Makefile. openmpi”  needs  to  be  modify.  Here  are  the  first  half  of  the  makefile  in  latest 
(Feb2010)  version  of  LAMMPS.  It  is  necessary  to  add  “CCFLAG”  and  “LINKFLAG” . 

jj  openmpi  =  Fedora  Core  6,  mpic++,  OpenMPI-1.1,  FFTW2 
SHELL  =  /bin/sh 

t - 

jj  compiler /linker  settings 

jj  specify  flags  and  libraries  needed  for  your  compiler 
CC  =  mpic++ 
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CCFLAGS  =  - 02  -funroll-loops  -DFFTJ'FTW  -DLAMMPS.GZIP 

-/strict- aliasing  -  Wall  -W  -W no-uninitialized 

-g  -0  -I/home/ (username) /include  -I/opt/openmpi.gcc/ 
include  -1/ opt/ openmpi.  gcc / include / openmpi 

DEP  FLAGS  =  -M 

LINK  =  mpic++ 

LINKFLAGS  =  -g  -0  -L/home/ (username) /lib  -L/opt/ 
openmpi. gcc /lib  -L/ opt/ openmpi. gcc/bin/ 

LIB  =  -Istdc-hP 

ARCHIVE  =  ar 

ARFLAGS  =  -rcsv 

SIZE  —  size 


This  file  is  updated  frequently  and  formats  are  slightly  different  in  each  version  of 
LAMMPS,  but  the  modifications  are  same. 

The  OpenMPI  is  used  for  parallelization  to  run  LAMMPS  and  the  command  is 
executed  in  a  shell  script. 

An  example  of  the  shell  script  is  as  following: 


cd  ( folder  name ) 

time  mpirun  ( execution  file )  <  ( input  file ) 

Here,  the  execution  file  for  OpenMPI  in  LAMMPS  is  “Imp. openmpi” .  The  shell 
script  file  needs  to  have  an  extension  of  “.sh”  and  it  is  run  by  “qsub”  command 
assigning  nodes  to  parallelize.  The  example  of  the  command  is: 
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qsub  —  l  nodes  =  (node  name )  :  ppn  =  ( number  of  CPU) 

+  ( node  name )  :  ppn  =  ( number  of  CPU)  ( shell  script  file) 

The  details  of  other  files  (input  file,  data  file  etc...)  can  be  found  in  the  LAMMPS 
manual  on  the  website  (http://lammps.sandia.gov/). 
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Appendix  C 


PIC  model 


PIC  simulation  has  been  done  to  try  to  observe  interactions  between  ion  jets.  The 
detail  descriptions  of  PIC  is  explained  in  M.S.  thesis  by  Fife  [81]. 

C.l  Code  Modification 

The  original  code  was  developed  by  Prof.  Paulo  Lozano  and  modifications  were 
applied  to  calculate  ion  emissions.  The  modified  code  was  able  to  increase  the  number 
of  sources  and  also  can  read  RPA  data  obtained  by  an  experiment.  The  RPA  data 
provided  an  information  about  the  fractions  of  solvated  ions  in  ion  jets.  We  use  RPA 
data  of  an  ionic  liquid  EMI-I  (l-Ethyl-3-methylimidazolium  Iodide)  obtained  by  Tim 
Fedkiw,  MIT  (figure  C-l). 

C.2  Two  Ion  Jets 

Figure  C-2  shows  visualization  of  the  simulation  results  in  chronological  order.  The 
simulation  time  step  is  5  [ns]  and  currents  for  monomer  and  dimer  are  63  [nA]  and  97 
[nA]  respectively.  The  beams  are  both  in  negative  polarity.  Interaction  can  be  seen 
from  jet  formations  which  are  apparently  not  straight  jets.  It  is  necessary  to  verify 
the  current  model  and  further  investigation  are  expected  to  adapt  the  model  to  array 
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Retarding  Potential  (-V) 


Figure  C-l:  RPA  data  for  EMI-I.  Courtesy  of  T.  Fedkiw  (MIT) 

of  electrosprays. 
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Figure  C-2:  Equilibrium  state  of  EMI-BF4  (a)  300K  (b)  350K  (c)  400K 
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